src/Core/TriangularMatrix.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
0005 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@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_TRIANGULARMATRIX_H
0012 #define EIGEN_TRIANGULARMATRIX_H
0013
0014 namespace Eigen {
0015
0016 namespace internal {
0017
0018 template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval;
0019
0020 }
0021
0022 /** \class TriangularBase
0023   * \ingroup Core_Module
0024   *
0025   * \brief Base class for triangular part in a matrix
0026   */
0027 template<typename Derived> class TriangularBase : public EigenBase<Derived>
0028 {
0029   public:
0030
0031     enum {
0032       Mode = internal::traits<Derived>::Mode,
0033       RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
0034       ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
0035       MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
0036       MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime,
0037
0038       SizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::RowsAtCompileTime,
0039                                                    internal::traits<Derived>::ColsAtCompileTime>::ret),
0040       /**< This is equal to the number of coefficients, i.e. the number of
0041           * rows times the number of columns, or to \a Dynamic if this is not
0042           * known at compile-time. \sa RowsAtCompileTime, ColsAtCompileTime */
0043
0044       MaxSizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::MaxRowsAtCompileTime,
0045                                                    internal::traits<Derived>::MaxColsAtCompileTime>::ret)
0046
0047     };
0048     typedef typename internal::traits<Derived>::Scalar Scalar;
0049     typedef typename internal::traits<Derived>::StorageKind StorageKind;
0050     typedef typename internal::traits<Derived>::StorageIndex StorageIndex;
0051     typedef typename internal::traits<Derived>::FullMatrixType DenseMatrixType;
0052     typedef DenseMatrixType DenseType;
0053     typedef Derived const& Nested;
0054
0055     EIGEN_DEVICE_FUNC
0056     inline TriangularBase() { eigen_assert(!((int(Mode) & int(UnitDiag)) && (int(Mode) & int(ZeroDiag)))); }
0057
0058     EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0059     inline Index rows() const EIGEN_NOEXCEPT { return derived().rows(); }
0060     EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0061     inline Index cols() const EIGEN_NOEXCEPT { return derived().cols(); }
0062     EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0063     inline Index outerStride() const EIGEN_NOEXCEPT { return derived().outerStride(); }
0064     EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0065     inline Index innerStride() const EIGEN_NOEXCEPT { return derived().innerStride(); }
0066
0067     // dummy resize function
0068     EIGEN_DEVICE_FUNC
0069     void resize(Index rows, Index cols)
0070     {
0071       EIGEN_UNUSED_VARIABLE(rows);
0072       EIGEN_UNUSED_VARIABLE(cols);
0073       eigen_assert(rows==this->rows() && cols==this->cols());
0074     }
0075
0076     EIGEN_DEVICE_FUNC
0077     inline Scalar coeff(Index row, Index col) const  { return derived().coeff(row,col); }
0078     EIGEN_DEVICE_FUNC
0079     inline Scalar& coeffRef(Index row, Index col) { return derived().coeffRef(row,col); }
0080
0081     /** \see MatrixBase::copyCoeff(row,col)
0082       */
0083     template<typename Other>
0084     EIGEN_DEVICE_FUNC
0085     EIGEN_STRONG_INLINE void copyCoeff(Index row, Index col, Other& other)
0086     {
0087       derived().coeffRef(row, col) = other.coeff(row, col);
0088     }
0089
0090     EIGEN_DEVICE_FUNC
0091     inline Scalar operator()(Index row, Index col) const
0092     {
0093       check_coordinates(row, col);
0094       return coeff(row,col);
0095     }
0096     EIGEN_DEVICE_FUNC
0097     inline Scalar& operator()(Index row, Index col)
0098     {
0099       check_coordinates(row, col);
0100       return coeffRef(row,col);
0101     }
0102
0103     #ifndef EIGEN_PARSED_BY_DOXYGEN
0104     EIGEN_DEVICE_FUNC
0105     inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
0106     EIGEN_DEVICE_FUNC
0107     inline Derived& derived() { return *static_cast<Derived*>(this); }
0108     #endif // not EIGEN_PARSED_BY_DOXYGEN
0109
0110     template<typename DenseDerived>
0111     EIGEN_DEVICE_FUNC
0112     void evalTo(MatrixBase<DenseDerived> &other) const;
0113     template<typename DenseDerived>
0114     EIGEN_DEVICE_FUNC
0115     void evalToLazy(MatrixBase<DenseDerived> &other) const;
0116
0117     EIGEN_DEVICE_FUNC
0118     DenseMatrixType toDenseMatrix() const
0119     {
0120       DenseMatrixType res(rows(), cols());
0121       evalToLazy(res);
0122       return res;
0123     }
0124
0125   protected:
0126
0127     void check_coordinates(Index row, Index col) const
0128     {
0129       EIGEN_ONLY_USED_FOR_DEBUG(row);
0130       EIGEN_ONLY_USED_FOR_DEBUG(col);
0131       eigen_assert(col>=0 && col<cols() && row>=0 && row<rows());
0132       const int mode = int(Mode) & ~SelfAdjoint;
0133       EIGEN_ONLY_USED_FOR_DEBUG(mode);
0134       eigen_assert((mode==Upper && col>=row)
0135                 || (mode==Lower && col<=row)
0136                 || ((mode==StrictlyUpper || mode==UnitUpper) && col>row)
0137                 || ((mode==StrictlyLower || mode==UnitLower) && col<row));
0138     }
0139
0140     #ifdef EIGEN_INTERNAL_DEBUGGING
0141     void check_coordinates_internal(Index row, Index col) const
0142     {
0143       check_coordinates(row, col);
0144     }
0145     #else
0146     void check_coordinates_internal(Index , Index ) const {}
0147     #endif
0148
0149 };
0150
0151 /** \class TriangularView
0152   * \ingroup Core_Module
0153   *
0154   * \brief Expression of a triangular part in a matrix
0155   *
0156   * \param MatrixType the type of the object in which we are taking the triangular part
0157   * \param Mode the kind of triangular matrix expression to construct. Can be #Upper,
0158   *             #Lower, #UnitUpper, #UnitLower, #StrictlyUpper, or #StrictlyLower.
0159   *             This is in fact a bit field; it must have either #Upper or #Lower,
0160   *             and additionally it may have #UnitDiag or #ZeroDiag or neither.
0161   *
0162   * This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular
0163   * matrices one should speak of "trapezoid" parts. This class is the return type
0164   * of MatrixBase::triangularView() and SparseMatrixBase::triangularView(), and most of the time this is the only way it is used.
0165   *
0166   * \sa MatrixBase::triangularView()
0167   */
0168 namespace internal {
0169 template<typename MatrixType, unsigned int _Mode>
0170 struct traits<TriangularView<MatrixType, _Mode> > : traits<MatrixType>
0171 {
0172   typedef typename ref_selector<MatrixType>::non_const_type MatrixTypeNested;
0173   typedef typename remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef;
0174   typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
0175   typedef typename MatrixType::PlainObject FullMatrixType;
0176   typedef MatrixType ExpressionType;
0177   enum {
0178     Mode = _Mode,
0179     FlagsLvalueBit = is_lvalue<MatrixType>::value ? LvalueBit : 0,
0180     Flags = (MatrixTypeNestedCleaned::Flags & (HereditaryBits | FlagsLvalueBit) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)))
0181   };
0182 };
0183 }
0184
0185 template<typename _MatrixType, unsigned int _Mode, typename StorageKind> class TriangularViewImpl;
0186
0187 template<typename _MatrixType, unsigned int _Mode> class TriangularView
0188   : public TriangularViewImpl<_MatrixType, _Mode, typename internal::traits<_MatrixType>::StorageKind >
0189 {
0190   public:
0191
0192     typedef TriangularViewImpl<_MatrixType, _Mode, typename internal::traits<_MatrixType>::StorageKind > Base;
0193     typedef typename internal::traits<TriangularView>::Scalar Scalar;
0194     typedef _MatrixType MatrixType;
0195
0196   protected:
0197     typedef typename internal::traits<TriangularView>::MatrixTypeNested MatrixTypeNested;
0198     typedef typename internal::traits<TriangularView>::MatrixTypeNestedNonRef MatrixTypeNestedNonRef;
0199
0200     typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType;
0201     typedef TriangularView<typename internal::add_const<MatrixType>::type, _Mode> ConstTriangularView;
0202
0203   public:
0204
0205     typedef typename internal::traits<TriangularView>::StorageKind StorageKind;
0206     typedef typename internal::traits<TriangularView>::MatrixTypeNestedCleaned NestedExpression;
0207
0208     enum {
0209       Mode = _Mode,
0210       Flags = internal::traits<TriangularView>::Flags,
0211       TransposeMode = (Mode & Upper ? Lower : 0)
0212                     | (Mode & Lower ? Upper : 0)
0213                     | (Mode & (UnitDiag))
0214                     | (Mode & (ZeroDiag)),
0215       IsVectorAtCompileTime = false
0216     };
0217
0218     EIGEN_DEVICE_FUNC
0219     explicit inline TriangularView(MatrixType& matrix) : m_matrix(matrix)
0220     {}
0221
0222     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(TriangularView)
0223
0224     /** \copydoc EigenBase::rows() */
0225     EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0226     inline Index rows() const EIGEN_NOEXCEPT { return m_matrix.rows(); }
0227     /** \copydoc EigenBase::cols() */
0228     EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0229     inline Index cols() const EIGEN_NOEXCEPT { return m_matrix.cols(); }
0230
0231     /** \returns a const reference to the nested expression */
0232     EIGEN_DEVICE_FUNC
0233     const NestedExpression& nestedExpression() const { return m_matrix; }
0234
0235     /** \returns a reference to the nested expression */
0236     EIGEN_DEVICE_FUNC
0237     NestedExpression& nestedExpression() { return m_matrix; }
0238
0239     typedef TriangularView<const MatrixConjugateReturnType,Mode> ConjugateReturnType;
0240     /** \sa MatrixBase::conjugate() const */
0241     EIGEN_DEVICE_FUNC
0242     inline const ConjugateReturnType conjugate() const
0243     { return ConjugateReturnType(m_matrix.conjugate()); }
0244
0245     /** \returns an expression of the complex conjugate of \c *this if Cond==true,
0246      *           returns \c *this otherwise.
0247      */
0248     template<bool Cond>
0249     EIGEN_DEVICE_FUNC
0250     inline typename internal::conditional<Cond,ConjugateReturnType,ConstTriangularView>::type
0251     conjugateIf() const
0252     {
0253       typedef typename internal::conditional<Cond,ConjugateReturnType,ConstTriangularView>::type ReturnType;
0254       return ReturnType(m_matrix.template conjugateIf<Cond>());
0255     }
0256
0257     typedef TriangularView<const typename MatrixType::AdjointReturnType,TransposeMode> AdjointReturnType;
0258     /** \sa MatrixBase::adjoint() const */
0259     EIGEN_DEVICE_FUNC
0260     inline const AdjointReturnType adjoint() const
0261     { return AdjointReturnType(m_matrix.adjoint()); }
0262
0263     typedef TriangularView<typename MatrixType::TransposeReturnType,TransposeMode> TransposeReturnType;
0264      /** \sa MatrixBase::transpose() */
0265     EIGEN_DEVICE_FUNC
0266     inline TransposeReturnType transpose()
0267     {
0268       EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
0269       typename MatrixType::TransposeReturnType tmp(m_matrix);
0270       return TransposeReturnType(tmp);
0271     }
0272
0273     typedef TriangularView<const typename MatrixType::ConstTransposeReturnType,TransposeMode> ConstTransposeReturnType;
0274     /** \sa MatrixBase::transpose() const */
0275     EIGEN_DEVICE_FUNC
0276     inline const ConstTransposeReturnType transpose() const
0277     {
0278       return ConstTransposeReturnType(m_matrix.transpose());
0279     }
0280
0281     template<typename Other>
0282     EIGEN_DEVICE_FUNC
0283     inline const Solve<TriangularView, Other>
0284     solve(const MatrixBase<Other>& other) const
0285     { return Solve<TriangularView, Other>(*this, other.derived()); }
0286
0287   // workaround MSVC ICE
0288   #if EIGEN_COMP_MSVC
0289     template<int Side, typename Other>
0290     EIGEN_DEVICE_FUNC
0291     inline const internal::triangular_solve_retval<Side,TriangularView, Other>
0292     solve(const MatrixBase<Other>& other) const
0293     { return Base::template solve<Side>(other); }
0294   #else
0295     using Base::solve;
0296   #endif
0297
0298     /** \returns a selfadjoint view of the referenced triangular part which must be either \c #Upper or \c #Lower.
0299       *
0300       * This is a shortcut for \code this->nestedExpression().selfadjointView<(*this)::Mode>() \endcode
0301       * \sa MatrixBase::selfadjointView() */
0302     EIGEN_DEVICE_FUNC
0303     SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView()
0304     {
0305       EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
0306       return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
0307     }
0308
0309     /** This is the const version of selfadjointView() */
0310     EIGEN_DEVICE_FUNC
0311     const SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() const
0312     {
0313       EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
0314       return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
0315     }
0316
0317
0318     /** \returns the determinant of the triangular matrix
0319       * \sa MatrixBase::determinant() */
0320     EIGEN_DEVICE_FUNC
0321     Scalar determinant() const
0322     {
0323       if (Mode & UnitDiag)
0324         return 1;
0325       else if (Mode & ZeroDiag)
0326         return 0;
0327       else
0328         return m_matrix.diagonal().prod();
0329     }
0330
0331   protected:
0332
0333     MatrixTypeNested m_matrix;
0334 };
0335
0336 /** \ingroup Core_Module
0337   *
0338   * \brief Base class for a triangular part in a \b dense matrix
0339   *
0340   * This class is an abstract base class of class TriangularView, and objects of type TriangularViewImpl cannot be instantiated.
0341   * It extends class TriangularView with additional methods which available for dense expressions only.
0342   *
0343   * \sa class TriangularView, MatrixBase::triangularView()
0344   */
0345 template<typename _MatrixType, unsigned int _Mode> class TriangularViewImpl<_MatrixType,_Mode,Dense>
0346   : public TriangularBase<TriangularView<_MatrixType, _Mode> >
0347 {
0348   public:
0349
0350     typedef TriangularView<_MatrixType, _Mode> TriangularViewType;
0351     typedef TriangularBase<TriangularViewType> Base;
0352     typedef typename internal::traits<TriangularViewType>::Scalar Scalar;
0353
0354     typedef _MatrixType MatrixType;
0355     typedef typename MatrixType::PlainObject DenseMatrixType;
0356     typedef DenseMatrixType PlainObject;
0357
0358   public:
0359     using Base::evalToLazy;
0360     using Base::derived;
0361
0362     typedef typename internal::traits<TriangularViewType>::StorageKind StorageKind;
0363
0364     enum {
0365       Mode = _Mode,
0366       Flags = internal::traits<TriangularViewType>::Flags
0367     };
0368
0369     /** \returns the outer-stride of the underlying dense matrix
0370       * \sa DenseCoeffsBase::outerStride() */
0371     EIGEN_DEVICE_FUNC
0372     inline Index outerStride() const { return derived().nestedExpression().outerStride(); }
0373     /** \returns the inner-stride of the underlying dense matrix
0374       * \sa DenseCoeffsBase::innerStride() */
0375     EIGEN_DEVICE_FUNC
0376     inline Index innerStride() const { return derived().nestedExpression().innerStride(); }
0377
0378     /** \sa MatrixBase::operator+=() */
0379     template<typename Other>
0380     EIGEN_DEVICE_FUNC
0381     TriangularViewType&  operator+=(const DenseBase<Other>& other) {
0382       internal::call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op<Scalar,typename Other::Scalar>());
0383       return derived();
0384     }
0385     /** \sa MatrixBase::operator-=() */
0386     template<typename Other>
0387     EIGEN_DEVICE_FUNC
0388     TriangularViewType&  operator-=(const DenseBase<Other>& other) {
0389       internal::call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op<Scalar,typename Other::Scalar>());
0390       return derived();
0391     }
0392
0393     /** \sa MatrixBase::operator*=() */
0394     EIGEN_DEVICE_FUNC
0395     TriangularViewType&  operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() * other; }
0396     /** \sa DenseBase::operator/=() */
0397     EIGEN_DEVICE_FUNC
0398     TriangularViewType&  operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() / other; }
0399
0400     /** \sa MatrixBase::fill() */
0401     EIGEN_DEVICE_FUNC
0402     void fill(const Scalar& value) { setConstant(value); }
0403     /** \sa MatrixBase::setConstant() */
0404     EIGEN_DEVICE_FUNC
0405     TriangularViewType& setConstant(const Scalar& value)
0406     { return *this = MatrixType::Constant(derived().rows(), derived().cols(), value); }
0407     /** \sa MatrixBase::setZero() */
0408     EIGEN_DEVICE_FUNC
0409     TriangularViewType& setZero() { return setConstant(Scalar(0)); }
0410     /** \sa MatrixBase::setOnes() */
0411     EIGEN_DEVICE_FUNC
0412     TriangularViewType& setOnes() { return setConstant(Scalar(1)); }
0413
0414     /** \sa MatrixBase::coeff()
0415       * \warning the coordinates must fit into the referenced triangular part
0416       */
0417     EIGEN_DEVICE_FUNC
0418     inline Scalar coeff(Index row, Index col) const
0419     {
0420       Base::check_coordinates_internal(row, col);
0421       return derived().nestedExpression().coeff(row, col);
0422     }
0423
0424     /** \sa MatrixBase::coeffRef()
0425       * \warning the coordinates must fit into the referenced triangular part
0426       */
0427     EIGEN_DEVICE_FUNC
0428     inline Scalar& coeffRef(Index row, Index col)
0429     {
0430       EIGEN_STATIC_ASSERT_LVALUE(TriangularViewType);
0431       Base::check_coordinates_internal(row, col);
0432       return derived().nestedExpression().coeffRef(row, col);
0433     }
0434
0435     /** Assigns a triangular matrix to a triangular part of a dense matrix */
0436     template<typename OtherDerived>
0437     EIGEN_DEVICE_FUNC
0438     TriangularViewType& operator=(const TriangularBase<OtherDerived>& other);
0439
0440     /** Shortcut for\code *this = other.other.triangularView<(*this)::Mode>() \endcode */
0441     template<typename OtherDerived>
0442     EIGEN_DEVICE_FUNC
0443     TriangularViewType& operator=(const MatrixBase<OtherDerived>& other);
0444
0445 #ifndef EIGEN_PARSED_BY_DOXYGEN
0446     EIGEN_DEVICE_FUNC
0447     TriangularViewType& operator=(const TriangularViewImpl& other)
0448     { return *this = other.derived().nestedExpression(); }
0449
0450     template<typename OtherDerived>
0451     /** \deprecated */
0452     EIGEN_DEPRECATED EIGEN_DEVICE_FUNC
0453     void lazyAssign(const TriangularBase<OtherDerived>& other);
0454
0455     template<typename OtherDerived>
0456     /** \deprecated */
0457     EIGEN_DEPRECATED EIGEN_DEVICE_FUNC
0458     void lazyAssign(const MatrixBase<OtherDerived>& other);
0459 #endif
0460
0461     /** Efficient triangular matrix times vector/matrix product */
0462     template<typename OtherDerived>
0463     EIGEN_DEVICE_FUNC
0464     const Product<TriangularViewType,OtherDerived>
0465     operator*(const MatrixBase<OtherDerived>& rhs) const
0466     {
0467       return Product<TriangularViewType,OtherDerived>(derived(), rhs.derived());
0468     }
0469
0470     /** Efficient vector/matrix times triangular matrix product */
0471     template<typename OtherDerived> friend
0472     EIGEN_DEVICE_FUNC
0473     const Product<OtherDerived,TriangularViewType>
0474     operator*(const MatrixBase<OtherDerived>& lhs, const TriangularViewImpl& rhs)
0475     {
0476       return Product<OtherDerived,TriangularViewType>(lhs.derived(),rhs.derived());
0477     }
0478
0479     /** \returns the product of the inverse of \c *this with \a other, \a *this being triangular.
0480       *
0481       * This function computes the inverse-matrix matrix product inverse(\c *this) * \a other if
0482       * \a Side==OnTheLeft (the default), or the right-inverse-multiply  \a other * inverse(\c *this) if
0483       * \a Side==OnTheRight.
0484       *
0485       * Note that the template parameter \c Side can be omitted, in which case \c Side==OnTheLeft
0486       *
0487       * The matrix \c *this must be triangular and invertible (i.e., all the coefficients of the
0488       * diagonal must be non zero). It works as a forward (resp. backward) substitution if \c *this
0489       * is an upper (resp. lower) triangular matrix.
0490       *
0491       * Example: \include Triangular_solve.cpp
0492       * Output: \verbinclude Triangular_solve.out
0493       *
0494       * This function returns an expression of the inverse-multiply and can works in-place if it is assigned
0495       * to the same matrix or vector \a other.
0496       *
0497       * For users coming from BLAS, this function (and more specifically solveInPlace()) offer
0498       * all the operations supported by the \c *TRSV and \c *TRSM BLAS routines.
0499       *
0500       * \sa TriangularView::solveInPlace()
0501       */
0502     template<int Side, typename Other>
0503     inline const internal::triangular_solve_retval<Side,TriangularViewType, Other>
0504     solve(const MatrixBase<Other>& other) const;
0505
0506     /** "in-place" version of TriangularView::solve() where the result is written in \a other
0507       *
0508       * \warning The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here.
0509       * This function will const_cast it, so constness isn't honored here.
0510       *
0511       * Note that the template parameter \c Side can be omitted, in which case \c Side==OnTheLeft
0512       *
0513       * See TriangularView:solve() for the details.
0514       */
0515     template<int Side, typename OtherDerived>
0516     EIGEN_DEVICE_FUNC
0517     void solveInPlace(const MatrixBase<OtherDerived>& other) const;
0518
0519     template<typename OtherDerived>
0520     EIGEN_DEVICE_FUNC
0521     void solveInPlace(const MatrixBase<OtherDerived>& other) const
0522     { return solveInPlace<OnTheLeft>(other); }
0523
0524     /** Swaps the coefficients of the common triangular parts of two matrices */
0525     template<typename OtherDerived>
0526     EIGEN_DEVICE_FUNC
0527 #ifdef EIGEN_PARSED_BY_DOXYGEN
0528     void swap(TriangularBase<OtherDerived> &other)
0529 #else
0530     void swap(TriangularBase<OtherDerived> const & other)
0531 #endif
0532     {
0533       EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
0534       call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
0535     }
0536
0537     /** Shortcut for \code (*this).swap(other.triangularView<(*this)::Mode>()) \endcode */
0538     template<typename OtherDerived>
0539     /** \deprecated */
0540     EIGEN_DEPRECATED EIGEN_DEVICE_FUNC
0541     void swap(MatrixBase<OtherDerived> const & other)
0542     {
0543       EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
0544       call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
0545     }
0546
0547     template<typename RhsType, typename DstType>
0548     EIGEN_DEVICE_FUNC
0549     EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const {
0550       if(!internal::is_same_dense(dst,rhs))
0551         dst = rhs;
0552       this->solveInPlace(dst);
0553     }
0554
0555     template<typename ProductType>
0556     EIGEN_DEVICE_FUNC
0557     EIGEN_STRONG_INLINE TriangularViewType& _assignProduct(const ProductType& prod, const Scalar& alpha, bool beta);
0558   protected:
0559     EIGEN_DEFAULT_COPY_CONSTRUCTOR(TriangularViewImpl)
0560     EIGEN_DEFAULT_EMPTY_CONSTRUCTOR_AND_DESTRUCTOR(TriangularViewImpl)
0561
0562 };
0563
0564 /***************************************************************************
0565 * Implementation of triangular evaluation/assignment
0566 ***************************************************************************/
0567
0568 #ifndef EIGEN_PARSED_BY_DOXYGEN
0569 // FIXME should we keep that possibility
0570 template<typename MatrixType, unsigned int Mode>
0571 template<typename OtherDerived>
0572 EIGEN_DEVICE_FUNC inline TriangularView<MatrixType, Mode>&
0573 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const MatrixBase<OtherDerived>& other)
0574 {
0575   internal::call_assignment_no_alias(derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
0576   return derived();
0577 }
0578
0579 // FIXME should we keep that possibility
0580 template<typename MatrixType, unsigned int Mode>
0581 template<typename OtherDerived>
0582 EIGEN_DEVICE_FUNC void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const MatrixBase<OtherDerived>& other)
0583 {
0584   internal::call_assignment_no_alias(derived(), other.template triangularView<Mode>());
0585 }
0586
0587
0588
0589 template<typename MatrixType, unsigned int Mode>
0590 template<typename OtherDerived>
0591 EIGEN_DEVICE_FUNC inline TriangularView<MatrixType, Mode>&
0592 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const TriangularBase<OtherDerived>& other)
0593 {
0594   eigen_assert(Mode == int(OtherDerived::Mode));
0595   internal::call_assignment(derived(), other.derived());
0596   return derived();
0597 }
0598
0599 template<typename MatrixType, unsigned int Mode>
0600 template<typename OtherDerived>
0601 EIGEN_DEVICE_FUNC void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const TriangularBase<OtherDerived>& other)
0602 {
0603   eigen_assert(Mode == int(OtherDerived::Mode));
0604   internal::call_assignment_no_alias(derived(), other.derived());
0605 }
0606 #endif
0607
0608 /***************************************************************************
0609 * Implementation of TriangularBase methods
0610 ***************************************************************************/
0611
0612 /** Assigns a triangular or selfadjoint matrix to a dense matrix.
0613   * If the matrix is triangular, the opposite part is set to zero. */
0614 template<typename Derived>
0615 template<typename DenseDerived>
0616 EIGEN_DEVICE_FUNC void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
0617 {
0618   evalToLazy(other.derived());
0619 }
0620
0621 /***************************************************************************
0622 * Implementation of TriangularView methods
0623 ***************************************************************************/
0624
0625 /***************************************************************************
0626 * Implementation of MatrixBase methods
0627 ***************************************************************************/
0628
0629 /**
0630   * \returns an expression of a triangular view extracted from the current matrix
0631   *
0632   * The parameter \a Mode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper,
0633   * \c #Lower, \c #StrictlyLower, \c #UnitLower.
0634   *
0635   * Example: \include MatrixBase_triangularView.cpp
0636   * Output: \verbinclude MatrixBase_triangularView.out
0637   *
0638   * \sa class TriangularView
0639   */
0640 template<typename Derived>
0641 template<unsigned int Mode>
0642 EIGEN_DEVICE_FUNC
0643 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
0644 MatrixBase<Derived>::triangularView()
0645 {
0646   return typename TriangularViewReturnType<Mode>::Type(derived());
0647 }
0648
0649 /** This is the const version of MatrixBase::triangularView() */
0650 template<typename Derived>
0651 template<unsigned int Mode>
0652 EIGEN_DEVICE_FUNC
0653 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
0654 MatrixBase<Derived>::triangularView() const
0655 {
0656   return typename ConstTriangularViewReturnType<Mode>::Type(derived());
0657 }
0658
0659 /** \returns true if *this is approximately equal to an upper triangular matrix,
0660   *          within the precision given by \a prec.
0661   *
0662   * \sa isLowerTriangular()
0663   */
0664 template<typename Derived>
0665 bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const
0666 {
0667   RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
0668   for(Index j = 0; j < cols(); ++j)
0669   {
0670     Index maxi = numext::mini(j, rows()-1);
0671     for(Index i = 0; i <= maxi; ++i)
0672     {
0673       RealScalar absValue = numext::abs(coeff(i,j));
0674       if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
0675     }
0676   }
0677   RealScalar threshold = maxAbsOnUpperPart * prec;
0678   for(Index j = 0; j < cols(); ++j)
0679     for(Index i = j+1; i < rows(); ++i)
0680       if(numext::abs(coeff(i, j)) > threshold) return false;
0681   return true;
0682 }
0683
0684 /** \returns true if *this is approximately equal to a lower triangular matrix,
0685   *          within the precision given by \a prec.
0686   *
0687   * \sa isUpperTriangular()
0688   */
0689 template<typename Derived>
0690 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const
0691 {
0692   RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
0693   for(Index j = 0; j < cols(); ++j)
0694     for(Index i = j; i < rows(); ++i)
0695     {
0696       RealScalar absValue = numext::abs(coeff(i,j));
0697       if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
0698     }
0699   RealScalar threshold = maxAbsOnLowerPart * prec;
0700   for(Index j = 1; j < cols(); ++j)
0701   {
0702     Index maxi = numext::mini(j, rows()-1);
0703     for(Index i = 0; i < maxi; ++i)
0704       if(numext::abs(coeff(i, j)) > threshold) return false;
0705   }
0706   return true;
0707 }
0708
0709
0710 /***************************************************************************
0711 ****************************************************************************
0712 * Evaluators and Assignment of triangular expressions
0713 ***************************************************************************
0714 ***************************************************************************/
0715
0716 namespace internal {
0717
0718
0719 // TODO currently a triangular expression has the form TriangularView<.,.>
0720 //      in the future triangular-ness should be defined by the expression traits
0721 //      such that Transpose<TriangularView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
0722 template<typename MatrixType, unsigned int Mode>
0723 struct evaluator_traits<TriangularView<MatrixType,Mode> >
0724 {
0725   typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
0726   typedef typename glue_shapes<typename evaluator_traits<MatrixType>::Shape, TriangularShape>::type Shape;
0727 };
0728
0729 template<typename MatrixType, unsigned int Mode>
0730 struct unary_evaluator<TriangularView<MatrixType,Mode>, IndexBased>
0731  : evaluator<typename internal::remove_all<MatrixType>::type>
0732 {
0733   typedef TriangularView<MatrixType,Mode> XprType;
0734   typedef evaluator<typename internal::remove_all<MatrixType>::type> Base;
0735   EIGEN_DEVICE_FUNC
0736   unary_evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {}
0737 };
0738
0739 // Additional assignment kinds:
0740 struct Triangular2Triangular    {};
0741 struct Triangular2Dense         {};
0742 struct Dense2Triangular         {};
0743
0744
0745 template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop;
0746
0747
0748 /** \internal Specialization of the dense assignment kernel for triangular matrices.
0749   * The main difference is that the triangular, diagonal, and opposite parts are processed through three different functions.
0750   * \tparam UpLo must be either Lower or Upper
0751   * \tparam Mode must be either 0, UnitDiag, ZeroDiag, or SelfAdjoint
0752   */
0753 template<int UpLo, int Mode, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version = Specialized>
0754 class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version>
0755 {
0756 protected:
0757   typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base;
0758   typedef typename Base::DstXprType DstXprType;
0759   typedef typename Base::SrcXprType SrcXprType;
0760   using Base::m_dst;
0761   using Base::m_src;
0762   using Base::m_functor;
0763 public:
0764
0765   typedef typename Base::DstEvaluatorType DstEvaluatorType;
0766   typedef typename Base::SrcEvaluatorType SrcEvaluatorType;
0767   typedef typename Base::Scalar Scalar;
0768   typedef typename Base::AssignmentTraits AssignmentTraits;
0769
0770
0771   EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr)
0772     : Base(dst, src, func, dstExpr)
0773   {}
0774
0775 #ifdef EIGEN_INTERNAL_DEBUGGING
0776   EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col)
0777   {
0778     eigen_internal_assert(row!=col);
0779     Base::assignCoeff(row,col);
0780   }
0781 #else
0782   using Base::assignCoeff;
0783 #endif
0784
0785   EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id)
0786   {
0787          if(Mode==UnitDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(1));
0788     else if(Mode==ZeroDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(0));
0789     else if(Mode==0)                       Base::assignCoeff(id,id);
0790   }
0791
0792   EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index row, Index col)
0793   {
0794     eigen_internal_assert(row!=col);
0795     if(SetOpposite)
0796       m_functor.assignCoeff(m_dst.coeffRef(row,col), Scalar(0));
0797   }
0798 };
0799
0800 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType, typename Functor>
0801 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0802 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src, const Functor &func)
0803 {
0804   typedef evaluator<DstXprType> DstEvaluatorType;
0805   typedef evaluator<SrcXprType> SrcEvaluatorType;
0806
0807   SrcEvaluatorType srcEvaluator(src);
0808
0809   Index dstRows = src.rows();
0810   Index dstCols = src.cols();
0811   if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
0812     dst.resize(dstRows, dstCols);
0813   DstEvaluatorType dstEvaluator(dst);
0814
0815   typedef triangular_dense_assignment_kernel< Mode&(Lower|Upper),Mode&(UnitDiag|ZeroDiag|SelfAdjoint),SetOpposite,
0816                                               DstEvaluatorType,SrcEvaluatorType,Functor> Kernel;
0817   Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived());
0818
0819   enum {
0820       unroll = DstXprType::SizeAtCompileTime != Dynamic
0821             && SrcEvaluatorType::CoeffReadCost < HugeCost
0822             && DstXprType::SizeAtCompileTime * (int(DstEvaluatorType::CoeffReadCost) + int(SrcEvaluatorType::CoeffReadCost)) / 2 <= EIGEN_UNROLLING_LIMIT
0823     };
0824
0825   triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, SetOpposite>::run(kernel);
0826 }
0827
0828 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType>
0829 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0830 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src)
0831 {
0832   call_triangular_assignment_loop<Mode,SetOpposite>(dst, src, internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>());
0833 }
0834
0835 template<> struct AssignmentKind<TriangularShape,TriangularShape> { typedef Triangular2Triangular Kind; };
0836 template<> struct AssignmentKind<DenseShape,TriangularShape>      { typedef Triangular2Dense      Kind; };
0837 template<> struct AssignmentKind<TriangularShape,DenseShape>      { typedef Dense2Triangular      Kind; };
0838
0839
0840 template< typename DstXprType, typename SrcXprType, typename Functor>
0841 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Triangular>
0842 {
0843   EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
0844   {
0845     eigen_assert(int(DstXprType::Mode) == int(SrcXprType::Mode));
0846
0847     call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
0848   }
0849 };
0850
0851 template< typename DstXprType, typename SrcXprType, typename Functor>
0852 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Dense>
0853 {
0854   EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
0855   {
0856     call_triangular_assignment_loop<SrcXprType::Mode, (int(SrcXprType::Mode) & int(SelfAdjoint)) == 0>(dst, src, func);
0857   }
0858 };
0859
0860 template< typename DstXprType, typename SrcXprType, typename Functor>
0861 struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular>
0862 {
0863   EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
0864   {
0865     call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
0866   }
0867 };
0868
0869
0870 template<typename Kernel, unsigned int Mode, int UnrollCount, bool SetOpposite>
0871 struct triangular_assignment_loop
0872 {
0873   // FIXME: this is not very clean, perhaps this information should be provided by the kernel?
0874   typedef typename Kernel::DstEvaluatorType DstEvaluatorType;
0875   typedef typename DstEvaluatorType::XprType DstXprType;
0876
0877   enum {
0878     col = (UnrollCount-1) / DstXprType::RowsAtCompileTime,
0879     row = (UnrollCount-1) % DstXprType::RowsAtCompileTime
0880   };
0881
0882   typedef typename Kernel::Scalar Scalar;
0883
0884   EIGEN_DEVICE_FUNC
0885   static inline void run(Kernel &kernel)
0886   {
0887     triangular_assignment_loop<Kernel, Mode, UnrollCount-1, SetOpposite>::run(kernel);
0888
0889     if(row==col)
0890       kernel.assignDiagonalCoeff(row);
0891     else if( ((Mode&Lower) && row>col) || ((Mode&Upper) && row<col) )
0892       kernel.assignCoeff(row,col);
0893     else if(SetOpposite)
0894       kernel.assignOppositeCoeff(row,col);
0895   }
0896 };
0897
0898 // prevent buggy user code from causing an infinite recursion
0899 template<typename Kernel, unsigned int Mode, bool SetOpposite>
0900 struct triangular_assignment_loop<Kernel, Mode, 0, SetOpposite>
0901 {
0902   EIGEN_DEVICE_FUNC
0903   static inline void run(Kernel &) {}
0904 };
0905
0906
0907
0908 // TODO: experiment with a recursive assignment procedure splitting the current
0909 //       triangular part into one rectangular and two triangular parts.
0910
0911
0912 template<typename Kernel, unsigned int Mode, bool SetOpposite>
0913 struct triangular_assignment_loop<Kernel, Mode, Dynamic, SetOpposite>
0914 {
0915   typedef typename Kernel::Scalar Scalar;
0916   EIGEN_DEVICE_FUNC
0917   static inline void run(Kernel &kernel)
0918   {
0919     for(Index j = 0; j < kernel.cols(); ++j)
0920     {
0921       Index maxi = numext::mini(j, kernel.rows());
0922       Index i = 0;
0923       if (((Mode&Lower) && SetOpposite) || (Mode&Upper))
0924       {
0925         for(; i < maxi; ++i)
0926           if(Mode&Upper) kernel.assignCoeff(i, j);
0927           else           kernel.assignOppositeCoeff(i, j);
0928       }
0929       else
0930         i = maxi;
0931
0932       if(i<kernel.rows()) // then i==j
0933         kernel.assignDiagonalCoeff(i++);
0934
0935       if (((Mode&Upper) && SetOpposite) || (Mode&Lower))
0936       {
0937         for(; i < kernel.rows(); ++i)
0938           if(Mode&Lower) kernel.assignCoeff(i, j);
0939           else           kernel.assignOppositeCoeff(i, j);
0940       }
0941     }
0942   }
0943 };
0944
0945 } // end namespace internal
0946
0947 /** Assigns a triangular or selfadjoint matrix to a dense matrix.
0948   * If the matrix is triangular, the opposite part is set to zero. */
0949 template<typename Derived>
0950 template<typename DenseDerived>
0951 EIGEN_DEVICE_FUNC void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
0952 {
0953   other.derived().resize(this->rows(), this->cols());
0954   internal::call_triangular_assignment_loop<Derived::Mode, (int(Derived::Mode) & int(SelfAdjoint)) == 0 /* SetOpposite */>(other.derived(), derived().nestedExpression());
0955 }
0956
0957 namespace internal {
0958
0959 // Triangular = Product
0960 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
0961 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
0962 {
0963   typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
0964   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,typename SrcXprType::Scalar> &)
0965   {
0966     Index dstRows = src.rows();
0967     Index dstCols = src.cols();
0968     if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
0969       dst.resize(dstRows, dstCols);
0970
0971     dst._assignProduct(src, Scalar(1), false);
0972   }
0973 };
0974
0975 // Triangular += Product
0976 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
0977 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::add_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
0978 {
0979   typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
0980   static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar,typename SrcXprType::Scalar> &)
0981   {
0982     dst._assignProduct(src, Scalar(1), true);
0983   }
0984 };
0985
0986 // Triangular -= Product
0987 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
0988 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::sub_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
0989 {
0990   typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
0991   static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,typename SrcXprType::Scalar> &)
0992   {
0993     dst._assignProduct(src, Scalar(-1), true);
0994   }
0995 };
0996
0997 } // end namespace internal
0998
0999 } // end namespace Eigen
1000
1001 #endif // EIGEN_TRIANGULARMATRIX_H