diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index fc6f56adc..09e960ae9 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -113,7 +113,7 @@ template class SparseSelfAdjointView SparseSelfAdjointView& rankUpdate(const SparseMatrixBase& u, Scalar alpha = Scalar(1)); /** \internal triggered by sparse_matrix = SparseSelfadjointView; */ - template void evalTo(SparseMatrix& _dest) const + template void evalTo(SparseMatrix& _dest) const { internal::permute_symm_to_fullsymm(m_matrix, _dest); } diff --git a/Eigen/src/SparseCore/SparseTriangularView.h b/Eigen/src/SparseCore/SparseTriangularView.h index 0b2d06528..3c0c10242 100644 --- a/Eigen/src/SparseCore/SparseTriangularView.h +++ b/Eigen/src/SparseCore/SparseTriangularView.h @@ -37,9 +37,10 @@ struct traits > template class SparseTriangularView : public SparseMatrixBase > { - enum { SkipFirst = (Mode==Lower && !(MatrixType::Flags&RowMajorBit)) - || (Mode==Upper && (MatrixType::Flags&RowMajorBit)), - SkipLast = !SkipFirst + enum { SkipFirst = ((Mode&Lower) && !(MatrixType::Flags&RowMajorBit)) + || ((Mode&Upper) && (MatrixType::Flags&RowMajorBit)), + SkipLast = !SkipFirst, + HasUnitDiag = (Mode&UnitDiag) ? 1 : 0 }; public: @@ -81,19 +82,61 @@ class SparseTriangularView::InnerIterator : public MatrixType:: public: EIGEN_STRONG_INLINE InnerIterator(const SparseTriangularView& view, Index outer) - : Base(view.nestedExpression(), outer) + : Base(view.nestedExpression(), outer), m_returnOne(false) { if(SkipFirst) - while((*this) && this->index()index()<=outer : this->index()=Base::outer())) + { + if((!SkipFirst) && Base::operator bool()) + Base::operator++(); + m_returnOne = true; + } } + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + if(HasUnitDiag && m_returnOne) + m_returnOne = false; + else + { + Base::operator++(); + if(HasUnitDiag && (!SkipFirst) && ((!Base::operator bool()) || Base::index()>=Base::outer())) + { + if((!SkipFirst) && Base::operator bool()) + Base::operator++(); + m_returnOne = true; + } + } + return *this; + } + inline Index row() const { return Base::row(); } inline Index col() const { return Base::col(); } + inline Index index() const + { + if(HasUnitDiag && m_returnOne) return Base::outer(); + else return Base::index(); + } + inline Scalar value() const + { + if(HasUnitDiag && m_returnOne) return Scalar(1); + else return Base::value(); + } EIGEN_STRONG_INLINE operator bool() const { - return SkipFirst ? Base::operator bool() : (Base::operator bool() && this->index() <= this->outer()); + if(HasUnitDiag && m_returnOne) + return true; + return (SkipFirst ? Base::operator bool() : (Base::operator bool() && this->index() <= this->outer())); } + protected: + bool m_returnOne; }; template @@ -105,10 +148,15 @@ class SparseTriangularView::ReverseInnerIterator : public Matri EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseTriangularView& view, Index outer) : Base(view.nestedExpression(), outer) { + eigen_assert((!HasUnitDiag) && "ReverseInnerIterator does not support yet triangular views with a unit diagonal"); if(SkipLast) while((*this) && this->index()>outer) --(*this); } + + EIGEN_STRONG_INLINE InnerIterator& operator--() + { Base::operator--(); return *this; } + inline Index row() const { return Base::row(); } inline Index col() const { return Base::col(); } diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index a3a6086cf..3a9cb61b3 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -196,7 +196,10 @@ template void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m1+=m2, refM1+=refM2); VERIFY_IS_APPROX(m1-=m2, refM1-=refM2); - VERIFY_IS_APPROX(m1.col(0).dot(refM2.row(0)), refM1.col(0).dot(refM2.row(0))); + if(SparseMatrixType::IsRowMajor) + VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.row(0).dot(refM2.row(0))); + else + VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.col(0).dot(refM2.row(0))); VERIFY_IS_APPROX(m1.conjugate(), refM1.conjugate()); VERIFY_IS_APPROX(m1.real(), refM1.real()); @@ -225,8 +228,15 @@ template void sparse_basic(const SparseMatrixType& re initSparse(density, refMat2, m2); int j0 = internal::random(0,rows-1); int j1 = internal::random(0,rows-1); - VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0)); - VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1)); + if(SparseMatrixType::IsRowMajor) + VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.row(j0)); + else + VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0)); + + if(SparseMatrixType::IsRowMajor) + VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.row(j0)+refMat2.row(j1)); + else + VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1)); //m2.innerVector(j0) = 2*m2.innerVector(j1); //refMat2.col(j0) = 2*refMat2.col(j1); //VERIFY_IS_APPROX(m2, refMat2); @@ -240,9 +250,16 @@ template void sparse_basic(const SparseMatrixType& re int j0 = internal::random(0,rows-2); int j1 = internal::random(0,rows-2); int n0 = internal::random(1,rows-(std::max)(j0,j1)); - VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0)); - VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0), - refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0)); + if(SparseMatrixType::IsRowMajor) + VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(j0,0,n0,cols)); + else + VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0)); + if(SparseMatrixType::IsRowMajor) + VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0), + refMat2.block(j0,0,n0,cols)+refMat2.block(j1,0,n0,cols)); + else + VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0), + refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0)); //m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0); //refMat2.block(0,j0,rows,n0) = refMat2.block(0,j0,rows,n0) + refMat2.block(0,j1,rows,n0); } @@ -272,7 +289,11 @@ template void sparse_basic(const SparseMatrixType& re else { countTrueNonZero++; - m2.insertBackByOuterInner(j,i) = refM2(i,j) = Scalar(1); + m2.insertBackByOuterInner(j,i) = Scalar(1); + if(SparseMatrixType::IsRowMajor) + refM2(j,i) = Scalar(1); + else + refM2(i,j) = Scalar(1); } } } @@ -283,8 +304,31 @@ template void sparse_basic(const SparseMatrixType& re VERIFY(countTrueNonZero==m2.nonZeros()); VERIFY_IS_APPROX(m2, refM2); } + + // test triangularView + { + DenseMatrix refMat2(rows, rows), refMat3(rows, rows); + SparseMatrixType m2(rows, rows), m3(rows, rows); + initSparse(density, refMat2, m2); + refMat3 = refMat2.template triangularView(); + m3 = m2.template triangularView(); + VERIFY_IS_APPROX(m3, refMat3); + + refMat3 = refMat2.template triangularView(); + m3 = m2.template triangularView(); + VERIFY_IS_APPROX(m3, refMat3); + + refMat3 = refMat2.template triangularView(); + m3 = m2.template triangularView(); + VERIFY_IS_APPROX(m3, refMat3); + + refMat3 = refMat2.template triangularView(); + m3 = m2.template triangularView(); + VERIFY_IS_APPROX(m3, refMat3); + } // test selfadjointView + if(!SparseMatrixType::IsRowMajor) { DenseMatrix refMat2(rows, rows), refMat3(rows, rows); SparseMatrixType m2(rows, rows), m3(rows, rows); @@ -308,8 +352,10 @@ void test_sparse_basic() for(int i = 0; i < g_repeat; i++) { int s = Eigen::internal::random(1,50); CALL_SUBTEST_1(( sparse_basic(SparseMatrix(8, 8)) )); - CALL_SUBTEST_2(( sparse_basic(SparseMatrix >(s, s)) )); + CALL_SUBTEST_2(( sparse_basic(SparseMatrix, ColMajor>(s, s)) )); + CALL_SUBTEST_2(( sparse_basic(SparseMatrix, RowMajor>(s, s)) )); CALL_SUBTEST_1(( sparse_basic(SparseMatrix(s, s)) )); CALL_SUBTEST_1(( sparse_basic(SparseMatrix(s, s)) )); + CALL_SUBTEST_1(( sparse_basic(SparseMatrix(s, s)) )); } }