diff options
Diffstat (limited to 'Eigen/src/SparseCore/SparseTriangularView.h')
| -rw-r--r-- | Eigen/src/SparseCore/SparseTriangularView.h | 252 |
1 files changed, 131 insertions, 121 deletions
diff --git a/Eigen/src/SparseCore/SparseTriangularView.h b/Eigen/src/SparseCore/SparseTriangularView.h index 333127b78..9ac120266 100644 --- a/Eigen/src/SparseCore/SparseTriangularView.h +++ b/Eigen/src/SparseCore/SparseTriangularView.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2009-2015 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla @@ -11,19 +11,19 @@ #ifndef EIGEN_SPARSE_TRIANGULARVIEW_H #define EIGEN_SPARSE_TRIANGULARVIEW_H -namespace Eigen { - -namespace internal { - -template<typename MatrixType, int Mode> -struct traits<SparseTriangularView<MatrixType,Mode> > -: public traits<MatrixType> -{}; - -} // namespace internal - -template<typename MatrixType, int Mode> class SparseTriangularView - : public SparseMatrixBase<SparseTriangularView<MatrixType,Mode> > +namespace Eigen { + +/** \ingroup SparseCore_Module + * + * \brief Base class for a triangular part in a \b sparse matrix + * + * This class is an abstract base class of class TriangularView, and objects of type TriangularViewImpl cannot be instantiated. + * It extends class TriangularView with additional methods which are available for sparse expressions only. + * + * \sa class TriangularView, SparseMatrixBase::triangularView() + */ +template<typename MatrixType, unsigned int Mode> class TriangularViewImpl<MatrixType,Mode,Sparse> + : public SparseMatrixBase<TriangularView<MatrixType,Mode> > { enum { SkipFirst = ((Mode&Lower) && !(MatrixType::Flags&RowMajorBit)) || ((Mode&Upper) && (MatrixType::Flags&RowMajorBit)), @@ -31,147 +31,157 @@ template<typename MatrixType, int Mode> class SparseTriangularView SkipDiag = (Mode&ZeroDiag) ? 1 : 0, HasUnitDiag = (Mode&UnitDiag) ? 1 : 0 }; + + typedef TriangularView<MatrixType,Mode> TriangularViewType; + + protected: + // dummy solve function to make TriangularView happy. + void solve() const; + typedef SparseMatrixBase<TriangularViewType> Base; 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(); } - + EIGEN_SPARSE_PUBLIC_INTERFACE(TriangularViewType) + typedef typename MatrixType::Nested MatrixTypeNested; typedef typename internal::remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef; typedef typename internal::remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned; - inline SparseTriangularView(const MatrixType& matrix) : m_matrix(matrix) {} - - /** \internal */ - inline const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; } - - template<typename OtherDerived> - typename internal::plain_matrix_type_column_major<OtherDerived>::type - solve(const MatrixBase<OtherDerived>& other) const; + template<typename RhsType, typename DstType> + EIGEN_DEVICE_FUNC + EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const { + if(!(internal::is_same<RhsType,DstType>::value && internal::extract_data(dst) == internal::extract_data(rhs))) + dst = rhs; + this->solveInPlace(dst); + } + /** Applies the inverse of \c *this to the dense vector or matrix \a other, "in-place" */ template<typename OtherDerived> void solveInPlace(MatrixBase<OtherDerived>& other) const; - template<typename OtherDerived> void solveInPlace(SparseMatrixBase<OtherDerived>& other) const; - protected: - MatrixTypeNested m_matrix; + /** Applies the inverse of \c *this to the sparse vector or matrix \a other, "in-place" */ + template<typename OtherDerived> void solveInPlace(SparseMatrixBase<OtherDerived>& other) const; + }; -template<typename MatrixType, int Mode> -class SparseTriangularView<MatrixType,Mode>::InnerIterator : public MatrixTypeNestedCleaned::InnerIterator +namespace internal { + +template<typename ArgType, unsigned int Mode> +struct unary_evaluator<TriangularView<ArgType,Mode>, IteratorBased> + : evaluator_base<TriangularView<ArgType,Mode> > { - typedef typename MatrixTypeNestedCleaned::InnerIterator Base; - typedef typename SparseTriangularView::Index Index; - public: + typedef TriangularView<ArgType,Mode> XprType; + +protected: + + typedef typename XprType::Scalar Scalar; + typedef typename XprType::StorageIndex StorageIndex; + typedef typename evaluator<ArgType>::InnerIterator EvalIterator; + + enum { SkipFirst = ((Mode&Lower) && !(ArgType::Flags&RowMajorBit)) + || ((Mode&Upper) && (ArgType::Flags&RowMajorBit)), + SkipLast = !SkipFirst, + SkipDiag = (Mode&ZeroDiag) ? 1 : 0, + HasUnitDiag = (Mode&UnitDiag) ? 1 : 0 + }; + +public: + + enum { + CoeffReadCost = evaluator<ArgType>::CoeffReadCost, + Flags = XprType::Flags + }; + + explicit unary_evaluator(const XprType &xpr) : m_argImpl(xpr.nestedExpression()), m_arg(xpr.nestedExpression()) {} + + inline Index nonZerosEstimate() const { + return m_argImpl.nonZerosEstimate(); + } + + class InnerIterator : public EvalIterator + { + typedef EvalIterator Base; + public: - EIGEN_STRONG_INLINE InnerIterator(const SparseTriangularView& view, Index outer) - : Base(view.nestedExpression(), outer), m_returnOne(false) - { - if(SkipFirst) + EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& xprEval, Index outer) + : Base(xprEval.m_argImpl,outer), m_returnOne(false), m_containsDiag(Base::outer()<xprEval.m_arg.innerSize()) { - while((*this) && ((HasUnitDiag||SkipDiag) ? this->index()<=outer : this->index()<outer)) - Base::operator++(); - if(HasUnitDiag) - m_returnOne = true; + if(SkipFirst) + { + while((*this) && ((HasUnitDiag||SkipDiag) ? this->index()<=outer : this->index()<outer)) + Base::operator++(); + if(HasUnitDiag) + m_returnOne = m_containsDiag; + } + else if(HasUnitDiag && ((!Base::operator bool()) || Base::index()>=Base::outer())) + { + if((!SkipFirst) && Base::operator bool()) + Base::operator++(); + m_returnOne = m_containsDiag; + } } - else if(HasUnitDiag && ((!Base::operator bool()) || Base::index()>=Base::outer())) + + EIGEN_STRONG_INLINE InnerIterator& operator++() { - if((!SkipFirst) && Base::operator bool()) + if(HasUnitDiag && m_returnOne) + m_returnOne = false; + else + { Base::operator++(); - m_returnOne = true; + if(HasUnitDiag && (!SkipFirst) && ((!Base::operator bool()) || Base::index()>=Base::outer())) + { + if((!SkipFirst) && Base::operator bool()) + Base::operator++(); + m_returnOne = m_containsDiag; + } + } + return *this; } - } - - EIGEN_STRONG_INLINE InnerIterator& operator++() - { - if(HasUnitDiag && m_returnOne) - m_returnOne = false; - else + + EIGEN_STRONG_INLINE operator bool() const { - Base::operator++(); - if(HasUnitDiag && (!SkipFirst) && ((!Base::operator bool()) || Base::index()>=Base::outer())) + if(HasUnitDiag && m_returnOne) + return true; + if(SkipFirst) return Base::operator bool(); + else { - if((!SkipFirst) && Base::operator bool()) - Base::operator++(); - m_returnOne = true; + if (SkipDiag) return (Base::operator bool() && this->index() < this->outer()); + else return (Base::operator bool() && this->index() <= this->outer()); } } - return *this; - } - - inline Index row() const { return (MatrixType::Flags&RowMajorBit ? Base::outer() : this->index()); } - inline Index col() const { return (MatrixType::Flags&RowMajorBit ? this->index() : Base::outer()); } - 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 - { - if(HasUnitDiag && m_returnOne) - return true; - if(SkipFirst) return Base::operator bool(); - else +// inline Index row() const { return (ArgType::Flags&RowMajorBit ? Base::outer() : this->index()); } +// inline Index col() const { return (ArgType::Flags&RowMajorBit ? this->index() : Base::outer()); } + inline StorageIndex index() const { - if (SkipDiag) return (Base::operator bool() && this->index() < this->outer()); - else return (Base::operator bool() && this->index() <= this->outer()); + if(HasUnitDiag && m_returnOne) return internal::convert_index<StorageIndex>(Base::outer()); + else return Base::index(); } - } - protected: - bool m_returnOne; -}; - -template<typename MatrixType, int Mode> -class SparseTriangularView<MatrixType,Mode>::ReverseInnerIterator : public MatrixTypeNestedCleaned::ReverseInnerIterator -{ - typedef typename MatrixTypeNestedCleaned::ReverseInnerIterator Base; - typedef typename SparseTriangularView::Index Index; - public: - - 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) && (SkipDiag ? this->index()>=outer : this->index()>outer)) - --(*this); - } - } - - EIGEN_STRONG_INLINE ReverseInnerIterator& operator--() - { Base::operator--(); return *this; } - - inline Index row() const { return Base::row(); } - inline Index col() const { return Base::col(); } - - EIGEN_STRONG_INLINE operator bool() const - { - if (SkipLast) return Base::operator bool() ; - else + inline Scalar value() const { - if(SkipDiag) return (Base::operator bool() && this->index() > this->outer()); - else return (Base::operator bool() && this->index() >= this->outer()); + if(HasUnitDiag && m_returnOne) return Scalar(1); + else return Base::value(); } - } + + protected: + bool m_returnOne; + bool m_containsDiag; + private: + Scalar& valueRef(); + }; + +protected: + evaluator<ArgType> m_argImpl; + const ArgType& m_arg; }; +} // end namespace internal + template<typename Derived> template<int Mode> -inline const SparseTriangularView<Derived, Mode> +inline const TriangularView<const Derived, Mode> SparseMatrixBase<Derived>::triangularView() const { - return derived(); + return TriangularView<const Derived, Mode>(derived()); } } // end namespace Eigen |