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