Main MRPT website > C++ reference
MRPT logo

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



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:16:28 UTC 2011