src/SparseCore/SparseTriangularView.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_SPARSE_TRIANGULARVIEW_H
0012 #define EIGEN_SPARSE_TRIANGULARVIEW_H
0013
0014 namespace Eigen {
0015
0016 /** \ingroup SparseCore_Module
0017   *
0018   * \brief Base class for a triangular part in a \b sparse matrix
0019   *
0020   * This class is an abstract base class of class TriangularView, and objects of type TriangularViewImpl cannot be instantiated.
0021   * It extends class TriangularView with additional methods which are available for sparse expressions only.
0022   *
0023   * \sa class TriangularView, SparseMatrixBase::triangularView()
0024   */
0025 template<typename MatrixType, unsigned int Mode> class TriangularViewImpl<MatrixType,Mode,Sparse>
0026   : public SparseMatrixBase<TriangularView<MatrixType,Mode> >
0027 {
0028     enum { SkipFirst = ((Mode&Lower) && !(MatrixType::Flags&RowMajorBit))
0029                     || ((Mode&Upper) &&  (MatrixType::Flags&RowMajorBit)),
0030            SkipLast = !SkipFirst,
0031            SkipDiag = (Mode&ZeroDiag) ? 1 : 0,
0032            HasUnitDiag = (Mode&UnitDiag) ? 1 : 0
0033     };
0034
0035     typedef TriangularView<MatrixType,Mode> TriangularViewType;
0036
0037   protected:
0038     // dummy solve function to make TriangularView happy.
0039     void solve() const;
0040
0041     typedef SparseMatrixBase<TriangularViewType> Base;
0042   public:
0043
0044     EIGEN_SPARSE_PUBLIC_INTERFACE(TriangularViewType)
0045
0046     typedef typename MatrixType::Nested MatrixTypeNested;
0047     typedef typename internal::remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef;
0048     typedef typename internal::remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
0049
0050     template<typename RhsType, typename DstType>
0051     EIGEN_DEVICE_FUNC
0052     EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const {
0053       if(!(internal::is_same<RhsType,DstType>::value && internal::extract_data(dst) == internal::extract_data(rhs)))
0054         dst = rhs;
0055       this->solveInPlace(dst);
0056     }
0057
0058     /** Applies the inverse of \c *this to the dense vector or matrix \a other, "in-place" */
0059     template<typename OtherDerived> void solveInPlace(MatrixBase<OtherDerived>& other) const;
0060
0061     /** Applies the inverse of \c *this to the sparse vector or matrix \a other, "in-place" */
0062     template<typename OtherDerived> void solveInPlace(SparseMatrixBase<OtherDerived>& other) const;
0063
0064 };
0065
0066 namespace internal {
0067
0068 template<typename ArgType, unsigned int Mode>
0069 struct unary_evaluator<TriangularView<ArgType,Mode>, IteratorBased>
0070  : evaluator_base<TriangularView<ArgType,Mode> >
0071 {
0072   typedef TriangularView<ArgType,Mode> XprType;
0073
0074 protected:
0075
0076   typedef typename XprType::Scalar Scalar;
0077   typedef typename XprType::StorageIndex StorageIndex;
0078   typedef typename evaluator<ArgType>::InnerIterator EvalIterator;
0079
0080   enum { SkipFirst = ((Mode&Lower) && !(ArgType::Flags&RowMajorBit))
0081                     || ((Mode&Upper) &&  (ArgType::Flags&RowMajorBit)),
0082          SkipLast = !SkipFirst,
0083          SkipDiag = (Mode&ZeroDiag) ? 1 : 0,
0084          HasUnitDiag = (Mode&UnitDiag) ? 1 : 0
0085   };
0086
0087 public:
0088
0089   enum {
0090     CoeffReadCost = evaluator<ArgType>::CoeffReadCost,
0091     Flags = XprType::Flags
0092   };
0093
0094   explicit unary_evaluator(const XprType &xpr) : m_argImpl(xpr.nestedExpression()), m_arg(xpr.nestedExpression()) {}
0095
0096   inline Index nonZerosEstimate() const {
0097     return m_argImpl.nonZerosEstimate();
0098   }
0099
0100   class InnerIterator : public EvalIterator
0101   {
0102       typedef EvalIterator Base;
0103     public:
0104
0105       EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& xprEval, Index outer)
0106         : Base(xprEval.m_argImpl,outer), m_returnOne(false), m_containsDiag(Base::outer()<xprEval.m_arg.innerSize())
0107       {
0108         if(SkipFirst)
0109         {
0110           while((*this) && ((HasUnitDiag||SkipDiag)  ? this->index()<=outer : this->index()<outer))
0111             Base::operator++();
0112           if(HasUnitDiag)
0113             m_returnOne = m_containsDiag;
0114         }
0115         else if(HasUnitDiag && ((!Base::operator bool()) || Base::index()>=Base::outer()))
0116         {
0117           if((!SkipFirst) && Base::operator bool())
0118             Base::operator++();
0119           m_returnOne = m_containsDiag;
0120         }
0121       }
0122
0123       EIGEN_STRONG_INLINE InnerIterator& operator++()
0124       {
0125         if(HasUnitDiag && m_returnOne)
0126           m_returnOne = false;
0127         else
0128         {
0129           Base::operator++();
0130           if(HasUnitDiag && (!SkipFirst) && ((!Base::operator bool()) || Base::index()>=Base::outer()))
0131           {
0132             if((!SkipFirst) && Base::operator bool())
0133               Base::operator++();
0134             m_returnOne = m_containsDiag;
0135           }
0136         }
0137         return *this;
0138       }
0139
0140       EIGEN_STRONG_INLINE operator bool() const
0141       {
0142         if(HasUnitDiag && m_returnOne)
0143           return true;
0144         if(SkipFirst) return  Base::operator bool();
0145         else
0146         {
0147           if (SkipDiag) return (Base::operator bool() && this->index() < this->outer());
0148           else return (Base::operator bool() && this->index() <= this->outer());
0149         }
0150       }
0151
0152 //       inline Index row() const { return (ArgType::Flags&RowMajorBit ? Base::outer() : this->index()); }
0153 //       inline Index col() const { return (ArgType::Flags&RowMajorBit ? this->index() : Base::outer()); }
0154       inline StorageIndex index() const
0155       {
0156         if(HasUnitDiag && m_returnOne)  return internal::convert_index<StorageIndex>(Base::outer());
0157         else                            return Base::index();
0158       }
0159       inline Scalar value() const
0160       {
0161         if(HasUnitDiag && m_returnOne)  return Scalar(1);
0162         else                            return Base::value();
0163       }
0164
0165     protected:
0166       bool m_returnOne;
0167       bool m_containsDiag;
0168     private:
0169       Scalar& valueRef();
0170   };
0171
0172 protected:
0173   evaluator<ArgType> m_argImpl;
0174   const ArgType& m_arg;
0175 };
0176
0177 } // end namespace internal
0178
0179 template<typename Derived>
0180 template<int Mode>
0181 inline const TriangularView<const Derived, Mode>
0182 SparseMatrixBase<Derived>::triangularView() const
0183 {
0184   return TriangularView<const Derived, Mode>(derived());
0185 }
0186
0187 } // end namespace Eigen
0188
0189 #endif // EIGEN_SPARSE_TRIANGULARVIEW_H