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


re_vision
Author(s): Dorian Galvez-Lopez
autogenerated on Sun Jan 5 2014 11:33:17