From 91e392a042ca8d40e460e5cf51d447bcce7a43d4 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 3 Dec 2011 23:49:37 +0100 Subject: add ReverseInnerIterators to loop over the elements in reverse order, and partly fix bug #356 (issue in trisolve for upper-column major)) --- Eigen/src/SparseCore/SparseCwiseUnaryOp.h | 93 +++++++++++++++++++++-------- Eigen/src/SparseCore/SparseMatrix.h | 34 +++++++++++ Eigen/src/SparseCore/SparseTriangularView.h | 28 ++++++++- Eigen/src/SparseCore/TriangularSolver.h | 8 ++- 4 files changed, 133 insertions(+), 30 deletions(-) (limited to 'Eigen/src/SparseCore') diff --git a/Eigen/src/SparseCore/SparseCwiseUnaryOp.h b/Eigen/src/SparseCore/SparseCwiseUnaryOp.h index a772c34bd..814649d38 100644 --- a/Eigen/src/SparseCore/SparseCwiseUnaryOp.h +++ b/Eigen/src/SparseCore/SparseCwiseUnaryOp.h @@ -32,39 +32,61 @@ class CwiseUnaryOpImpl public: class InnerIterator; -// typedef typename internal::remove_reference::type _LhsNested; + class ReverseInnerIterator; typedef CwiseUnaryOp Derived; EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) + + protected: + typedef typename internal::traits::_XprTypeNested _MatrixTypeNested; + typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator; + typedef typename _MatrixTypeNested::ReverseInnerIterator MatrixTypeReverseIterator; }; template class CwiseUnaryOpImpl::InnerIterator + : public CwiseUnaryOpImpl::MatrixTypeIterator { typedef typename CwiseUnaryOpImpl::Scalar Scalar; - typedef typename internal::traits::_XprTypeNested _MatrixTypeNested; - typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator; - typedef typename MatrixType::Index Index; + typedef typename CwiseUnaryOpImpl::MatrixTypeIterator Base; public: EIGEN_STRONG_INLINE InnerIterator(const CwiseUnaryOpImpl& unaryOp, Index outer) - : m_iter(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) + : Base(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) {} EIGEN_STRONG_INLINE InnerIterator& operator++() - { ++m_iter; return *this; } + { Base::operator++(); return *this; } + + EIGEN_STRONG_INLINE Scalar value() const { return m_functor(Base::value()); } - EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_iter.value()); } + protected: + const UnaryOp m_functor; + private: + Scalar& valueRef(); +}; + +template +class CwiseUnaryOpImpl::ReverseInnerIterator + : public CwiseUnaryOpImpl::MatrixTypeReverseIterator +{ + typedef typename CwiseUnaryOpImpl::Scalar Scalar; + typedef typename CwiseUnaryOpImpl::MatrixTypeReverseIterator Base; + public: + + EIGEN_STRONG_INLINE ReverseInnerIterator(const CwiseUnaryOpImpl& unaryOp, Index outer) + : Base(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) + {} - EIGEN_STRONG_INLINE Index index() const { return m_iter.index(); } - EIGEN_STRONG_INLINE Index row() const { return m_iter.row(); } - EIGEN_STRONG_INLINE Index col() const { return m_iter.col(); } + EIGEN_STRONG_INLINE ReverseInnerIterator& operator--() + { Base::operator--(); return *this; } - EIGEN_STRONG_INLINE operator bool() const { return m_iter; } + EIGEN_STRONG_INLINE Scalar value() const { return m_functor(Base::value()); } protected: - MatrixTypeIterator m_iter; const UnaryOp m_functor; + private: + Scalar& valueRef(); }; template @@ -74,39 +96,58 @@ class CwiseUnaryViewImpl public: class InnerIterator; -// typedef typename internal::remove_reference::type _LhsNested; + class ReverseInnerIterator; typedef CwiseUnaryView Derived; EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) + + protected: + typedef typename internal::traits::_MatrixTypeNested _MatrixTypeNested; + typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator; + typedef typename _MatrixTypeNested::ReverseInnerIterator MatrixTypeReverseIterator; }; template class CwiseUnaryViewImpl::InnerIterator + : public CwiseUnaryViewImpl::MatrixTypeIterator { typedef typename CwiseUnaryViewImpl::Scalar Scalar; - typedef typename internal::traits::_MatrixTypeNested _MatrixTypeNested; - typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator; - typedef typename MatrixType::Index Index; + typedef typename CwiseUnaryViewImpl::MatrixTypeIterator Base; public: - EIGEN_STRONG_INLINE InnerIterator(const CwiseUnaryViewImpl& unaryView, Index outer) - : m_iter(unaryView.derived().nestedExpression(),outer), m_functor(unaryView.derived().functor()) + EIGEN_STRONG_INLINE InnerIterator(const CwiseUnaryViewImpl& unaryOp, Index outer) + : Base(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) {} EIGEN_STRONG_INLINE InnerIterator& operator++() - { ++m_iter; return *this; } + { Base::operator++(); return *this; } + + EIGEN_STRONG_INLINE Scalar value() const { return m_functor(Base::value()); } + EIGEN_STRONG_INLINE Scalar& valueRef() { return m_functor(Base::valueRef()); } - EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_iter.value()); } - EIGEN_STRONG_INLINE Scalar& valueRef() { return m_functor(m_iter.valueRef()); } + protected: + const ViewOp m_functor; +}; + +template +class CwiseUnaryViewImpl::ReverseInnerIterator + : public CwiseUnaryViewImpl::MatrixTypeReverseIterator +{ + typedef typename CwiseUnaryViewImpl::Scalar Scalar; + typedef typename CwiseUnaryViewImpl::MatrixTypeReverseIterator Base; + public: + + EIGEN_STRONG_INLINE ReverseInnerIterator(const CwiseUnaryViewImpl& unaryOp, Index outer) + : Base(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) + {} - EIGEN_STRONG_INLINE Index index() const { return m_iter.index(); } - EIGEN_STRONG_INLINE Index row() const { return m_iter.row(); } - EIGEN_STRONG_INLINE Index col() const { return m_iter.col(); } + EIGEN_STRONG_INLINE ReverseInnerIterator& operator--() + { Base::operator--(); return *this; } - EIGEN_STRONG_INLINE operator bool() const { return m_iter; } + EIGEN_STRONG_INLINE Scalar value() const { return m_functor(Base::value()); } + EIGEN_STRONG_INLINE Scalar& valueRef() { return m_functor(Base::valueRef()); } protected: - MatrixTypeIterator m_iter; const ViewOp m_functor; }; diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 2d938a671..cc2805a53 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -210,6 +210,7 @@ class SparseMatrix public: class InnerIterator; + class ReverseInnerIterator; /** Removes all non zeros but keep allocated memory */ inline void setZero() @@ -889,4 +890,37 @@ class SparseMatrix::InnerIterator Index m_end; }; +template +class SparseMatrix::ReverseInnerIterator +{ + public: + ReverseInnerIterator(const SparseMatrix& mat, Index outer) + : m_values(mat._valuePtr()), m_indices(mat._innerIndexPtr()), m_outer(outer), m_start(mat.m_outerIndex[outer]) + { + if(mat.compressed()) + m_id = mat.m_outerIndex[outer+1]; + else + m_id = m_start + mat.m_innerNonZeros[outer]; + } + + inline ReverseInnerIterator& operator--() { --m_id; return *this; } + + inline const Scalar& value() const { return m_values[m_id-1]; } + inline Scalar& valueRef() { return const_cast(m_values[m_id-1]); } + + inline Index index() const { return m_indices[m_id-1]; } + inline Index outer() const { return m_outer; } + inline Index row() const { return IsRowMajor ? m_outer : index(); } + inline Index col() const { return IsRowMajor ? index() : m_outer; } + + inline operator bool() const { return (m_id > m_start); } + + protected: + const Scalar* m_values; + const Index* m_indices; + const Index m_outer; + Index m_id; + const Index m_start; +}; + #endif // EIGEN_SPARSEMATRIX_H diff --git a/Eigen/src/SparseCore/SparseTriangularView.h b/Eigen/src/SparseCore/SparseTriangularView.h index dff5ae2c4..0b2d06528 100644 --- a/Eigen/src/SparseCore/SparseTriangularView.h +++ b/Eigen/src/SparseCore/SparseTriangularView.h @@ -38,12 +38,16 @@ template class SparseTriangularView : public SparseMatrixBase > { enum { SkipFirst = (Mode==Lower && !(MatrixType::Flags&RowMajorBit)) - || (Mode==Upper && (MatrixType::Flags&RowMajorBit)) }; + || (Mode==Upper && (MatrixType::Flags&RowMajorBit)), + SkipLast = !SkipFirst + }; + public: EIGEN_SPARSE_PUBLIC_INTERFACE(SparseTriangularView) class InnerIterator; + class ReverseInnerIterator; inline Index rows() const { return m_matrix.rows(); } inline Index cols() const { return m_matrix.cols(); } @@ -92,6 +96,28 @@ class SparseTriangularView::InnerIterator : public MatrixType:: } }; +template +class SparseTriangularView::ReverseInnerIterator : public MatrixType::ReverseInnerIterator +{ + typedef typename MatrixType::ReverseInnerIterator Base; + public: + + EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseTriangularView& view, Index outer) + : Base(view.nestedExpression(), outer) + { + if(SkipLast) + while((*this) && this->index()>outer) + --(*this); + } + inline Index row() const { return Base::row(); } + inline Index col() const { return Base::col(); } + + EIGEN_STRONG_INLINE operator bool() const + { + return SkipLast ? Base::operator bool() : (Base::operator bool() && this->index() >= this->outer()); + } +}; + template template inline const SparseTriangularView diff --git a/Eigen/src/SparseCore/TriangularSolver.h b/Eigen/src/SparseCore/TriangularSolver.h index 6c66052bf..45128afe2 100644 --- a/Eigen/src/SparseCore/TriangularSolver.h +++ b/Eigen/src/SparseCore/TriangularSolver.h @@ -156,9 +156,11 @@ struct sparse_solve_triangular_selector { if(!(Mode & UnitDiag)) { - // FIXME lhs.coeff(i,i) might not be always efficient while it must simply be the - // last element of the column ! - other.coeffRef(i,col) /= lhs.innerVector(i).lastCoeff(); + typename Lhs::ReverseInnerIterator it(lhs, i); + while(it && it.index()!=i) + --it; + eigen_assert(it && it.index()==i); + other.coeffRef(i,col) /= it.value(); } typename Lhs::InnerIterator it(lhs, i); for(; it && it.index()