TriangularMatrix.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_TRIANGULARMATRIX_H
00012 #define EIGEN_TRIANGULARMATRIX_H
00013 
00014 namespace Eigen { 
00015 
00016 namespace internal {
00017   
00018 template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval;
00019   
00020 }
00021 
00022 /** \internal
00023   *
00024   * \class TriangularBase
00025   * \ingroup Core_Module
00026   *
00027   * \brief Base class for triangular part in a matrix
00028   */
00029 template<typename Derived> class TriangularBase : public EigenBase<Derived>
00030 {
00031   public:
00032 
00033     enum {
00034       Mode = internal::traits<Derived>::Mode,
00035       CoeffReadCost = internal::traits<Derived>::CoeffReadCost,
00036       RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
00037       ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
00038       MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
00039       MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime
00040     };
00041     typedef typename internal::traits<Derived>::Scalar Scalar;
00042     typedef typename internal::traits<Derived>::StorageKind StorageKind;
00043     typedef typename internal::traits<Derived>::Index Index;
00044     typedef typename internal::traits<Derived>::DenseMatrixType DenseMatrixType;
00045     typedef DenseMatrixType DenseType;
00046 
00047     inline TriangularBase() { eigen_assert(!((Mode&UnitDiag) && (Mode&ZeroDiag))); }
00048 
00049     inline Index rows() const { return derived().rows(); }
00050     inline Index cols() const { return derived().cols(); }
00051     inline Index outerStride() const { return derived().outerStride(); }
00052     inline Index innerStride() const { return derived().innerStride(); }
00053 
00054     inline Scalar coeff(Index row, Index col) const  { return derived().coeff(row,col); }
00055     inline Scalar& coeffRef(Index row, Index col) { return derived().coeffRef(row,col); }
00056 
00057     /** \see MatrixBase::copyCoeff(row,col)
00058       */
00059     template<typename Other>
00060     EIGEN_STRONG_INLINE void copyCoeff(Index row, Index col, Other& other)
00061     {
00062       derived().coeffRef(row, col) = other.coeff(row, col);
00063     }
00064 
00065     inline Scalar operator()(Index row, Index col) const
00066     {
00067       check_coordinates(row, col);
00068       return coeff(row,col);
00069     }
00070     inline Scalar& operator()(Index row, Index col)
00071     {
00072       check_coordinates(row, col);
00073       return coeffRef(row,col);
00074     }
00075 
00076     #ifndef EIGEN_PARSED_BY_DOXYGEN
00077     inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
00078     inline Derived& derived() { return *static_cast<Derived*>(this); }
00079     #endif // not EIGEN_PARSED_BY_DOXYGEN
00080 
00081     template<typename DenseDerived>
00082     void evalTo(MatrixBase<DenseDerived> &other) const;
00083     template<typename DenseDerived>
00084     void evalToLazy(MatrixBase<DenseDerived> &other) const;
00085 
00086     DenseMatrixType toDenseMatrix() const
00087     {
00088       DenseMatrixType res(rows(), cols());
00089       evalToLazy(res);
00090       return res;
00091     }
00092 
00093   protected:
00094 
00095     void check_coordinates(Index row, Index col) const
00096     {
00097       EIGEN_ONLY_USED_FOR_DEBUG(row);
00098       EIGEN_ONLY_USED_FOR_DEBUG(col);
00099       eigen_assert(col>=0 && col<cols() && row>=0 && row<rows());
00100       const int mode = int(Mode) & ~SelfAdjoint;
00101       EIGEN_ONLY_USED_FOR_DEBUG(mode);
00102       eigen_assert((mode==Upper && col>=row)
00103                 || (mode==Lower && col<=row)
00104                 || ((mode==StrictlyUpper || mode==UnitUpper) && col>row)
00105                 || ((mode==StrictlyLower || mode==UnitLower) && col<row));
00106     }
00107 
00108     #ifdef EIGEN_INTERNAL_DEBUGGING
00109     void check_coordinates_internal(Index row, Index col) const
00110     {
00111       check_coordinates(row, col);
00112     }
00113     #else
00114     void check_coordinates_internal(Index , Index ) const {}
00115     #endif
00116 
00117 };
00118 
00119 /** \class TriangularView
00120   * \ingroup Core_Module
00121   *
00122   * \brief Base class for triangular part in a matrix
00123   *
00124   * \param MatrixType the type of the object in which we are taking the triangular part
00125   * \param Mode the kind of triangular matrix expression to construct. Can be #Upper,
00126   *             #Lower, #UnitUpper, #UnitLower, #StrictlyUpper, or #StrictlyLower.
00127   *             This is in fact a bit field; it must have either #Upper or #Lower, 
00128   *             and additionnaly it may have #UnitDiag or #ZeroDiag or neither.
00129   *
00130   * This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular
00131   * matrices one should speak of "trapezoid" parts. This class is the return type
00132   * of MatrixBase::triangularView() and most of the time this is the only way it is used.
00133   *
00134   * \sa MatrixBase::triangularView()
00135   */
00136 namespace internal {
00137 template<typename MatrixType, unsigned int _Mode>
00138 struct traits<TriangularView<MatrixType, _Mode> > : traits<MatrixType>
00139 {
00140   typedef typename nested<MatrixType>::type MatrixTypeNested;
00141   typedef typename remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef;
00142   typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
00143   typedef MatrixType ExpressionType;
00144   typedef typename MatrixType::PlainObject DenseMatrixType;
00145   enum {
00146     Mode = _Mode,
00147     Flags = (MatrixTypeNestedCleaned::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) | Mode,
00148     CoeffReadCost = MatrixTypeNestedCleaned::CoeffReadCost
00149   };
00150 };
00151 }
00152 
00153 template<int Mode, bool LhsIsTriangular,
00154          typename Lhs, bool LhsIsVector,
00155          typename Rhs, bool RhsIsVector>
00156 struct TriangularProduct;
00157 
00158 template<typename _MatrixType, unsigned int _Mode> class TriangularView
00159   : public TriangularBase<TriangularView<_MatrixType, _Mode> >
00160 {
00161   public:
00162 
00163     typedef TriangularBase<TriangularView> Base;
00164     typedef typename internal::traits<TriangularView>::Scalar Scalar;
00165 
00166     typedef _MatrixType MatrixType;
00167     typedef typename internal::traits<TriangularView>::DenseMatrixType DenseMatrixType;
00168     typedef DenseMatrixType PlainObject;
00169 
00170   protected:
00171     typedef typename internal::traits<TriangularView>::MatrixTypeNested MatrixTypeNested;
00172     typedef typename internal::traits<TriangularView>::MatrixTypeNestedNonRef MatrixTypeNestedNonRef;
00173     typedef typename internal::traits<TriangularView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned;
00174 
00175     typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType;
00176     
00177   public:
00178     using Base::evalToLazy;
00179   
00180 
00181     typedef typename internal::traits<TriangularView>::StorageKind StorageKind;
00182     typedef typename internal::traits<TriangularView>::Index Index;
00183 
00184     enum {
00185       Mode = _Mode,
00186       TransposeMode = (Mode & Upper ? Lower : 0)
00187                     | (Mode & Lower ? Upper : 0)
00188                     | (Mode & (UnitDiag))
00189                     | (Mode & (ZeroDiag))
00190     };
00191 
00192     inline TriangularView(const MatrixType& matrix) : m_matrix(matrix)
00193     {}
00194 
00195     inline Index rows () const { return m_matrix.rows(); }
00196     inline Index cols () const { return m_matrix.cols(); }
00197     inline Index outerStride() const { return m_matrix.outerStride(); }
00198     inline Index innerStride() const { return m_matrix.innerStride(); }
00199 
00200     /** \sa MatrixBase::operator+=() */
00201     template<typename Other> TriangularView&  operator+= (const DenseBase<Other>& other) { return *this = m_matrix + other.derived(); }
00202     /** \sa MatrixBase::operator-=() */
00203     template<typename Other> TriangularView&  operator-= (const DenseBase<Other>& other) { return *this = m_matrix - other.derived(); }
00204     /** \sa MatrixBase::operator*=() */
00205     TriangularView&  operator*= (const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix * other; }
00206     /** \sa MatrixBase::operator/=() */
00207     TriangularView&  operator/= (const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix / other; }
00208 
00209     /** \sa MatrixBase::fill() */
00210     void fill (const Scalar& value) { setConstant (value); }
00211     /** \sa MatrixBase::setConstant() */
00212     TriangularView& setConstant (const Scalar& value)
00213     { return *this = MatrixType::Constant(rows (), cols (), value); }
00214     /** \sa MatrixBase::setZero() */
00215     TriangularView& setZero () { return setConstant (Scalar(0)); }
00216     /** \sa MatrixBase::setOnes() */
00217     TriangularView& setOnes () { return setConstant (Scalar(1)); }
00218 
00219     /** \sa MatrixBase::coeff()
00220       * \warning the coordinates must fit into the referenced triangular part
00221       */
00222     inline Scalar coeff (Index row, Index col) const
00223     {
00224       Base::check_coordinates_internal(row, col);
00225       return m_matrix.coeff(row, col);
00226     }
00227 
00228     /** \sa MatrixBase::coeffRef()
00229       * \warning the coordinates must fit into the referenced triangular part
00230       */
00231     inline Scalar& coeffRef (Index row, Index col)
00232     {
00233       Base::check_coordinates_internal(row, col);
00234       return m_matrix.const_cast_derived().coeffRef(row, col);
00235     }
00236 
00237     const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; }
00238     MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); }
00239 
00240     /** Assigns a triangular matrix to a triangular part of a dense matrix */
00241     template<typename OtherDerived>
00242     TriangularView& operator=(const TriangularBase<OtherDerived>& other);
00243 
00244     template<typename OtherDerived>
00245     TriangularView& operator=(const MatrixBase<OtherDerived>& other);
00246 
00247     TriangularView& operator=(const TriangularView& other)
00248     { return *this = other.nestedExpression(); }
00249 
00250     template<typename OtherDerived>
00251     void lazyAssign(const TriangularBase<OtherDerived>& other);
00252 
00253     template<typename OtherDerived>
00254     void lazyAssign(const MatrixBase<OtherDerived>& other);
00255 
00256     /** \sa MatrixBase::conjugate() */
00257     inline TriangularView<MatrixConjugateReturnType,Mode> conjugate ()
00258     { return m_matrix.conjugate(); }
00259     /** \sa MatrixBase::conjugate() const */
00260     inline const TriangularView<MatrixConjugateReturnType,Mode> conjugate () const
00261     { return m_matrix.conjugate(); }
00262 
00263     /** \sa MatrixBase::adjoint() const */
00264     inline const TriangularView<const typename MatrixType::AdjointReturnType,TransposeMode> adjoint () const
00265     { return m_matrix.adjoint(); }
00266 
00267     /** \sa MatrixBase::transpose() */
00268     inline TriangularView<Transpose<MatrixType>,TransposeMode> transpose ()
00269     {
00270       EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
00271       return m_matrix.const_cast_derived().transpose();
00272     }
00273     /** \sa MatrixBase::transpose() const */
00274     inline const TriangularView<Transpose<MatrixType>,TransposeMode> transpose () const
00275     {
00276       return m_matrix.transpose();
00277     }
00278 
00279     /** Efficient triangular matrix times vector/matrix product */
00280     template<typename OtherDerived>
00281     TriangularProduct<Mode, true, MatrixType, false, OtherDerived, OtherDerived::ColsAtCompileTime==1>
00282     operator*(const MatrixBase<OtherDerived>& rhs) const
00283     {
00284       return TriangularProduct
00285               <Mode, true, MatrixType, false, OtherDerived, OtherDerived::ColsAtCompileTime==1>
00286               (m_matrix, rhs.derived());
00287     }
00288 
00289     /** Efficient vector/matrix times triangular matrix product */
00290     template<typename OtherDerived> friend
00291     TriangularProduct<Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false>
00292     operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs)
00293     {
00294       return TriangularProduct
00295               <Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false>
00296               (lhs.derived(),rhs.m_matrix);
00297     }
00298 
00299     #ifdef EIGEN2_SUPPORT
00300     template<typename OtherDerived>
00301     struct eigen2_product_return_type
00302     {
00303       typedef typename TriangularView<MatrixType,Mode>::DenseMatrixType DenseMatrixType;
00304       typedef typename OtherDerived::PlainObject::DenseType OtherPlainObject;
00305       typedef typename ProductReturnType<DenseMatrixType, OtherPlainObject>::Type ProdRetType;
00306       typedef typename ProdRetType::PlainObject type;
00307     };
00308     template<typename OtherDerived>
00309     const typename eigen2_product_return_type<OtherDerived>::type
00310     operator*(const EigenBase<OtherDerived>& rhs) const
00311     {
00312       typename OtherDerived::PlainObject::DenseType rhsPlainObject;
00313       rhs.evalTo(rhsPlainObject);
00314       return this->toDenseMatrix() * rhsPlainObject;
00315     }
00316     template<typename OtherMatrixType>
00317     bool isApprox(const TriangularView<OtherMatrixType, Mode>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const
00318     {
00319       return this->toDenseMatrix().isApprox(other.toDenseMatrix(), precision);
00320     }
00321     template<typename OtherDerived>
00322     bool isApprox(const MatrixBase<OtherDerived>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const
00323     {
00324       return this->toDenseMatrix().isApprox(other, precision);
00325     }
00326     #endif // EIGEN2_SUPPORT
00327 
00328     template<int Side, typename Other>
00329     inline const internal::triangular_solve_retval<Side,TriangularView, Other>
00330     solve (const MatrixBase<Other>& other) const;
00331 
00332     template<int Side, typename OtherDerived>
00333     void solveInPlace(const MatrixBase<OtherDerived>& other) const;
00334 
00335     template<typename Other>
00336     inline const internal::triangular_solve_retval<OnTheLeft,TriangularView, Other> 
00337     solve (const MatrixBase<Other>& other) const
00338     { return solve<OnTheLeft>(other); }
00339 
00340     template<typename OtherDerived>
00341     void solveInPlace(const MatrixBase<OtherDerived>& other) const
00342     { return solveInPlace<OnTheLeft>(other); }
00343 
00344     const SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() const
00345     {
00346       EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
00347       return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
00348     }
00349     SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView()
00350     {
00351       EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
00352       return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
00353     }
00354 
00355     template<typename OtherDerived>
00356     void swap(TriangularBase<OtherDerived> const & other)
00357     {
00358       TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived());
00359     }
00360 
00361     template<typename OtherDerived>
00362     void swap(MatrixBase<OtherDerived> const & other)
00363     {
00364       SwapWrapper<MatrixType> swaper(const_cast<MatrixType&>(m_matrix));
00365       TriangularView<SwapWrapper<MatrixType>,Mode>(swaper).lazyAssign(other.derived());
00366     }
00367 
00368     Scalar determinant() const
00369     {
00370       if (Mode & UnitDiag)
00371         return 1;
00372       else if (Mode & ZeroDiag)
00373         return 0;
00374       else
00375         return m_matrix.diagonal().prod();
00376     }
00377     
00378     // TODO simplify the following:
00379     template<typename ProductDerived, typename Lhs, typename Rhs>
00380     EIGEN_STRONG_INLINE TriangularView& operator=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
00381     {
00382       setZero ();
00383       return assignProduct(other.derived(),1);
00384     }
00385     
00386     template<typename ProductDerived, typename Lhs, typename Rhs>
00387     EIGEN_STRONG_INLINE TriangularView& operator+= (const ProductBase<ProductDerived, Lhs,Rhs>& other)
00388     {
00389       return assignProduct(other.derived(),1);
00390     }
00391     
00392     template<typename ProductDerived, typename Lhs, typename Rhs>
00393     EIGEN_STRONG_INLINE TriangularView& operator-= (const ProductBase<ProductDerived, Lhs,Rhs>& other)
00394     {
00395       return assignProduct(other.derived(),-1);
00396     }
00397     
00398     
00399     template<typename ProductDerived>
00400     EIGEN_STRONG_INLINE TriangularView& operator=(const ScaledProduct<ProductDerived>& other)
00401     {
00402       setZero ();
00403       return assignProduct(other.derived(),other.alpha());
00404     }
00405     
00406     template<typename ProductDerived>
00407     EIGEN_STRONG_INLINE TriangularView& operator+= (const ScaledProduct<ProductDerived>& other)
00408     {
00409       return assignProduct(other.derived(),other.alpha());
00410     }
00411     
00412     template<typename ProductDerived>
00413     EIGEN_STRONG_INLINE TriangularView& operator-= (const ScaledProduct<ProductDerived>& other)
00414     {
00415       return assignProduct(other.derived(),-other.alpha());
00416     }
00417     
00418   protected:
00419     
00420     template<typename ProductDerived, typename Lhs, typename Rhs>
00421     EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase<ProductDerived, Lhs,Rhs>& prod, const Scalar& alpha);
00422     
00423     template<int Mode, bool LhsIsTriangular,
00424          typename Lhs, bool LhsIsVector,
00425          typename Rhs, bool RhsIsVector>
00426     EIGEN_STRONG_INLINE TriangularView& assignProduct(const TriangularProduct<Mode, LhsIsTriangular, Lhs, LhsIsVector, Rhs, RhsIsVector>& prod, const Scalar& alpha)
00427     {
00428       lazyAssign(alpha*prod.eval());
00429       return *this;
00430     }
00431 
00432     MatrixTypeNested m_matrix;
00433 };
00434 
00435 /***************************************************************************
00436 * Implementation of triangular evaluation/assignment
00437 ***************************************************************************/
00438 
00439 namespace internal {
00440 
00441 template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount, bool ClearOpposite>
00442 struct triangular_assignment_selector
00443 {
00444   enum {
00445     col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
00446     row = (UnrollCount-1) % Derived1::RowsAtCompileTime
00447   };
00448   
00449   typedef typename Derived1::Scalar Scalar;
00450 
00451   static inline void run(Derived1 &dst, const Derived2 &src)
00452   {
00453     triangular_assignment_selector<Derived1, Derived2, Mode, UnrollCount-1, ClearOpposite>::run(dst, src);
00454 
00455     eigen_assert( Mode == Upper || Mode == Lower
00456             || Mode == StrictlyUpper || Mode == StrictlyLower
00457             || Mode == UnitUpper || Mode == UnitLower);
00458     if((Mode == Upper && row <= col)
00459     || (Mode == Lower && row >= col)
00460     || (Mode == StrictlyUpper && row < col)
00461     || (Mode == StrictlyLower && row > col)
00462     || (Mode == UnitUpper && row < col)
00463     || (Mode == UnitLower && row > col))
00464       dst.copyCoeff(row, col, src);
00465     else if(ClearOpposite)
00466     {
00467       if (Mode&UnitDiag && row==col)
00468         dst.coeffRef(row, col) = Scalar(1);
00469       else
00470         dst.coeffRef(row, col) = Scalar(0);
00471     }
00472   }
00473 };
00474 
00475 // prevent buggy user code from causing an infinite recursion
00476 template<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite>
00477 struct triangular_assignment_selector<Derived1, Derived2, Mode, 0, ClearOpposite>
00478 {
00479   static inline void run(Derived1 &, const Derived2 &) {}
00480 };
00481 
00482 template<typename Derived1, typename Derived2, bool ClearOpposite>
00483 struct triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearOpposite>
00484 {
00485   typedef typename Derived1::Index Index;
00486   typedef typename Derived1::Scalar Scalar;
00487   static inline void run(Derived1 &dst, const Derived2 &src)
00488   {
00489     for(Index j = 0; j < dst.cols(); ++j)
00490     {
00491       Index maxi = (std::min)(j, dst.rows()-1);
00492       for(Index i = 0; i <= maxi; ++i)
00493         dst.copyCoeff(i, j, src);
00494       if (ClearOpposite)
00495         for(Index i = maxi+1; i < dst.rows(); ++i)
00496           dst.coeffRef(i, j) = Scalar(0);
00497     }
00498   }
00499 };
00500 
00501 template<typename Derived1, typename Derived2, bool ClearOpposite>
00502 struct triangular_assignment_selector<Derived1, Derived2, Lower, Dynamic, ClearOpposite>
00503 {
00504   typedef typename Derived1::Index Index;
00505   static inline void run(Derived1 &dst, const Derived2 &src)
00506   {
00507     for(Index j = 0; j < dst.cols(); ++j)
00508     {
00509       for(Index i = j; i < dst.rows(); ++i)
00510         dst.copyCoeff(i, j, src);
00511       Index maxi = (std::min)(j, dst.rows());
00512       if (ClearOpposite)
00513         for(Index i = 0; i < maxi; ++i)
00514           dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
00515     }
00516   }
00517 };
00518 
00519 template<typename Derived1, typename Derived2, bool ClearOpposite>
00520 struct triangular_assignment_selector<Derived1, Derived2, StrictlyUpper, Dynamic, ClearOpposite>
00521 {
00522   typedef typename Derived1::Index Index;
00523   typedef typename Derived1::Scalar Scalar;
00524   static inline void run(Derived1 &dst, const Derived2 &src)
00525   {
00526     for(Index j = 0; j < dst.cols(); ++j)
00527     {
00528       Index maxi = (std::min)(j, dst.rows());
00529       for(Index i = 0; i < maxi; ++i)
00530         dst.copyCoeff(i, j, src);
00531       if (ClearOpposite)
00532         for(Index i = maxi; i < dst.rows(); ++i)
00533           dst.coeffRef(i, j) = Scalar(0);
00534     }
00535   }
00536 };
00537 
00538 template<typename Derived1, typename Derived2, bool ClearOpposite>
00539 struct triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic, ClearOpposite>
00540 {
00541   typedef typename Derived1::Index Index;
00542   static inline void run(Derived1 &dst, const Derived2 &src)
00543   {
00544     for(Index j = 0; j < dst.cols(); ++j)
00545     {
00546       for(Index i = j+1; i < dst.rows(); ++i)
00547         dst.copyCoeff(i, j, src);
00548       Index maxi = (std::min)(j, dst.rows()-1);
00549       if (ClearOpposite)
00550         for(Index i = 0; i <= maxi; ++i)
00551           dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
00552     }
00553   }
00554 };
00555 
00556 template<typename Derived1, typename Derived2, bool ClearOpposite>
00557 struct triangular_assignment_selector<Derived1, Derived2, UnitUpper, Dynamic, ClearOpposite>
00558 {
00559   typedef typename Derived1::Index Index;
00560   static inline void run(Derived1 &dst, const Derived2 &src)
00561   {
00562     for(Index j = 0; j < dst.cols(); ++j)
00563     {
00564       Index maxi = (std::min)(j, dst.rows());
00565       for(Index i = 0; i < maxi; ++i)
00566         dst.copyCoeff(i, j, src);
00567       if (ClearOpposite)
00568       {
00569         for(Index i = maxi+1; i < dst.rows(); ++i)
00570           dst.coeffRef(i, j) = 0;
00571       }
00572     }
00573     dst.diagonal().setOnes();
00574   }
00575 };
00576 template<typename Derived1, typename Derived2, bool ClearOpposite>
00577 struct triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, ClearOpposite>
00578 {
00579   typedef typename Derived1::Index Index;
00580   static inline void run(Derived1 &dst, const Derived2 &src)
00581   {
00582     for(Index j = 0; j < dst.cols(); ++j)
00583     {
00584       Index maxi = (std::min)(j, dst.rows());
00585       for(Index i = maxi+1; i < dst.rows(); ++i)
00586         dst.copyCoeff(i, j, src);
00587       if (ClearOpposite)
00588       {
00589         for(Index i = 0; i < maxi; ++i)
00590           dst.coeffRef(i, j) = 0;
00591       }
00592     }
00593     dst.diagonal().setOnes();
00594   }
00595 };
00596 
00597 } // end namespace internal
00598 
00599 // FIXME should we keep that possibility
00600 template<typename MatrixType, unsigned int Mode>
00601 template<typename OtherDerived>
00602 inline TriangularView<MatrixType, Mode>&
00603 TriangularView<MatrixType, Mode>::operator=(const MatrixBase<OtherDerived>& other)
00604 {
00605   if(OtherDerived::Flags & EvalBeforeAssigningBit)
00606   {
00607     typename internal::plain_matrix_type<OtherDerived>::type other_evaluated(other.rows(), other.cols());
00608     other_evaluated.template triangularView<Mode>().lazyAssign(other.derived());
00609     lazyAssign(other_evaluated);
00610   }
00611   else
00612     lazyAssign(other.derived());
00613   return *this;
00614 }
00615 
00616 // FIXME should we keep that possibility
00617 template<typename MatrixType, unsigned int Mode>
00618 template<typename OtherDerived>
00619 void TriangularView<MatrixType, Mode>::lazyAssign(const MatrixBase<OtherDerived>& other)
00620 {
00621   enum {
00622     unroll = MatrixType::SizeAtCompileTime != Dynamic
00623           && internal::traits<OtherDerived>::CoeffReadCost != Dynamic
00624           && MatrixType::SizeAtCompileTime*internal::traits<OtherDerived>::CoeffReadCost/2 <= EIGEN_UNROLLING_LIMIT
00625   };
00626   eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
00627 
00628   internal::triangular_assignment_selector
00629     <MatrixType, OtherDerived, int(Mode),
00630     unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic,
00631     false // do not change the opposite triangular part
00632     >::run(m_matrix.const_cast_derived(), other.derived());
00633 }
00634 
00635 
00636 
00637 template<typename MatrixType, unsigned int Mode>
00638 template<typename OtherDerived>
00639 inline TriangularView<MatrixType, Mode>&
00640 TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>& other)
00641 {
00642   eigen_assert(Mode == int(OtherDerived::Mode));
00643   if(internal::traits<OtherDerived>::Flags & EvalBeforeAssigningBit)
00644   {
00645     typename OtherDerived::DenseMatrixType other_evaluated(other.rows(), other.cols());
00646     other_evaluated.template triangularView<Mode>().lazyAssign(other.derived().nestedExpression());
00647     lazyAssign(other_evaluated);
00648   }
00649   else
00650     lazyAssign(other.derived().nestedExpression());
00651   return *this;
00652 }
00653 
00654 template<typename MatrixType, unsigned int Mode>
00655 template<typename OtherDerived>
00656 void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDerived>& other)
00657 {
00658   enum {
00659     unroll = MatrixType::SizeAtCompileTime != Dynamic
00660                    && internal::traits<OtherDerived>::CoeffReadCost != Dynamic
00661                    && MatrixType::SizeAtCompileTime * internal::traits<OtherDerived>::CoeffReadCost / 2
00662                         <= EIGEN_UNROLLING_LIMIT
00663   };
00664   eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
00665 
00666   internal::triangular_assignment_selector
00667     <MatrixType, OtherDerived, int(Mode),
00668     unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic,
00669     false // preserve the opposite triangular part
00670     >::run(m_matrix.const_cast_derived(), other.derived().nestedExpression());
00671 }
00672 
00673 /***************************************************************************
00674 * Implementation of TriangularBase methods
00675 ***************************************************************************/
00676 
00677 /** Assigns a triangular or selfadjoint matrix to a dense matrix.
00678   * If the matrix is triangular, the opposite part is set to zero. */
00679 template<typename Derived>
00680 template<typename DenseDerived>
00681 void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
00682 {
00683   if(internal::traits<Derived>::Flags & EvalBeforeAssigningBit)
00684   {
00685     typename internal::plain_matrix_type<Derived>::type other_evaluated(rows(), cols());
00686     evalToLazy(other_evaluated);
00687     other.derived().swap(other_evaluated);
00688   }
00689   else
00690     evalToLazy(other.derived());
00691 }
00692 
00693 /** Assigns a triangular or selfadjoint matrix to a dense matrix.
00694   * If the matrix is triangular, the opposite part is set to zero. */
00695 template<typename Derived>
00696 template<typename DenseDerived>
00697 void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
00698 {
00699   enum {
00700     unroll = DenseDerived::SizeAtCompileTime != Dynamic
00701                    && internal::traits<Derived>::CoeffReadCost != Dynamic
00702                    && DenseDerived::SizeAtCompileTime * internal::traits<Derived>::CoeffReadCost / 2
00703                         <= EIGEN_UNROLLING_LIMIT
00704   };
00705   other.derived().resize(this->rows(), this->cols());
00706 
00707   internal::triangular_assignment_selector
00708     <DenseDerived, typename internal::traits<Derived>::MatrixTypeNestedCleaned, Derived::Mode,
00709     unroll ? int(DenseDerived::SizeAtCompileTime) : Dynamic,
00710     true // clear the opposite triangular part
00711     >::run(other.derived(), derived().nestedExpression());
00712 }
00713 
00714 /***************************************************************************
00715 * Implementation of TriangularView methods
00716 ***************************************************************************/
00717 
00718 /***************************************************************************
00719 * Implementation of MatrixBase methods
00720 ***************************************************************************/
00721 
00722 #ifdef EIGEN2_SUPPORT
00723 
00724 // implementation of part<>(), including the SelfAdjoint case.
00725 
00726 namespace internal {
00727 template<typename MatrixType, unsigned int Mode>
00728 struct eigen2_part_return_type
00729 {
00730   typedef TriangularView<MatrixType, Mode> type;
00731 };
00732 
00733 template<typename MatrixType>
00734 struct eigen2_part_return_type<MatrixType, SelfAdjoint>
00735 {
00736   typedef SelfAdjointView<MatrixType, Upper> type;
00737 };
00738 }
00739 
00740 /** \deprecated use MatrixBase::triangularView() */
00741 template<typename Derived>
00742 template<unsigned int Mode>
00743 const typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part () const
00744 {
00745   return derived();
00746 }
00747 
00748 /** \deprecated use MatrixBase::triangularView() */
00749 template<typename Derived>
00750 template<unsigned int Mode>
00751 typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part ()
00752 {
00753   return derived();
00754 }
00755 #endif
00756 
00757 /**
00758   * \returns an expression of a triangular view extracted from the current matrix
00759   *
00760   * The parameter \a Mode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper,
00761   * \c #Lower, \c #StrictlyLower, \c #UnitLower.
00762   *
00763   * Example: \include MatrixBase_extract.cpp
00764   * Output: \verbinclude MatrixBase_extract.out
00765   *
00766   * \sa class TriangularView
00767   */
00768 template<typename Derived>
00769 template<unsigned int Mode>
00770 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
00771 MatrixBase<Derived>::triangularView ()
00772 {
00773   return derived();
00774 }
00775 
00776 /** This is the const version of MatrixBase::triangularView() */
00777 template<typename Derived>
00778 template<unsigned int Mode>
00779 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
00780 MatrixBase<Derived>::triangularView () const
00781 {
00782   return derived();
00783 }
00784 
00785 /** \returns true if *this is approximately equal to an upper triangular matrix,
00786   *          within the precision given by \a prec.
00787   *
00788   * \sa isLowerTriangular()
00789   */
00790 template<typename Derived>
00791 bool MatrixBase<Derived>::isUpperTriangular (const RealScalar& prec) const
00792 {
00793   using std::abs;
00794   RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
00795   for(Index j = 0; j < cols(); ++j)
00796   {
00797     Index maxi = (std::min)(j, rows()-1);
00798     for(Index i = 0; i <= maxi; ++i)
00799     {
00800       RealScalar absValue = abs(coeff(i,j));
00801       if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
00802     }
00803   }
00804   RealScalar threshold = maxAbsOnUpperPart * prec;
00805   for(Index j = 0; j < cols(); ++j)
00806     for(Index i = j+1; i < rows(); ++i)
00807       if(abs(coeff(i, j)) > threshold) return false;
00808   return true;
00809 }
00810 
00811 /** \returns true if *this is approximately equal to a lower triangular matrix,
00812   *          within the precision given by \a prec.
00813   *
00814   * \sa isUpperTriangular()
00815   */
00816 template<typename Derived>
00817 bool MatrixBase<Derived>::isLowerTriangular (const RealScalar& prec) const
00818 {
00819   using std::abs;
00820   RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
00821   for(Index j = 0; j < cols(); ++j)
00822     for(Index i = j; i < rows(); ++i)
00823     {
00824       RealScalar absValue = abs(coeff(i,j));
00825       if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
00826     }
00827   RealScalar threshold = maxAbsOnLowerPart * prec;
00828   for(Index j = 1; j < cols(); ++j)
00829   {
00830     Index maxi = (std::min)(j, rows()-1);
00831     for(Index i = 0; i < maxi; ++i)
00832       if(abs(coeff(i, j)) > threshold) return false;
00833   }
00834   return true;
00835 }
00836 
00837 } // end namespace Eigen
00838 
00839 #endif // EIGEN_TRIANGULARMATRIX_H