Eigen  3.3.3
TriangularMatrix.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_TRIANGULARMATRIX_H
00012 #define EIGEN_TRIANGULARMATRIX_H
00013 
00014 namespace Eigen { 
00015 
00016 namespace internal {
00017   
00018 template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval;
00019   
00020 }
00021 
00027 template<typename Derived> class TriangularBase : public EigenBase<Derived>
00028 {
00029   public:
00030 
00031     enum {
00032       Mode = internal::traits<Derived>::Mode,
00033       RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
00034       ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
00035       MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
00036       MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime,
00037       
00038       SizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::RowsAtCompileTime,
00039                                                    internal::traits<Derived>::ColsAtCompileTime>::ret),
00044       MaxSizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::MaxRowsAtCompileTime,
00045                                                    internal::traits<Derived>::MaxColsAtCompileTime>::ret)
00046         
00047     };
00048     typedef typename internal::traits<Derived>::Scalar Scalar;
00049     typedef typename internal::traits<Derived>::StorageKind StorageKind;
00050     typedef typename internal::traits<Derived>::StorageIndex StorageIndex;
00051     typedef typename internal::traits<Derived>::FullMatrixType DenseMatrixType;
00052     typedef DenseMatrixType DenseType;
00053     typedef Derived const& Nested;
00054 
00055     EIGEN_DEVICE_FUNC
00056     inline TriangularBase() { eigen_assert(!((Mode&UnitDiag) && (Mode&ZeroDiag))); }
00057 
00058     EIGEN_DEVICE_FUNC
00059     inline Index rows() const { return derived().rows(); }
00060     EIGEN_DEVICE_FUNC
00061     inline Index cols() const { return derived().cols(); }
00062     EIGEN_DEVICE_FUNC
00063     inline Index outerStride() const { return derived().outerStride(); }
00064     EIGEN_DEVICE_FUNC
00065     inline Index innerStride() const { return derived().innerStride(); }
00066     
00067     // dummy resize function
00068     void resize(Index rows, Index cols)
00069     {
00070       EIGEN_UNUSED_VARIABLE(rows);
00071       EIGEN_UNUSED_VARIABLE(cols);
00072       eigen_assert(rows==this->rows() && cols==this->cols());
00073     }
00074 
00075     EIGEN_DEVICE_FUNC
00076     inline Scalar coeff(Index row, Index col) const  { return derived().coeff(row,col); }
00077     EIGEN_DEVICE_FUNC
00078     inline Scalar& coeffRef(Index row, Index col) { return derived().coeffRef(row,col); }
00079 
00082     template<typename Other>
00083     EIGEN_DEVICE_FUNC
00084     EIGEN_STRONG_INLINE void copyCoeff(Index row, Index col, Other& other)
00085     {
00086       derived().coeffRef(row, col) = other.coeff(row, col);
00087     }
00088 
00089     EIGEN_DEVICE_FUNC
00090     inline Scalar operator()(Index row, Index col) const
00091     {
00092       check_coordinates(row, col);
00093       return coeff(row,col);
00094     }
00095     EIGEN_DEVICE_FUNC
00096     inline Scalar& operator()(Index row, Index col)
00097     {
00098       check_coordinates(row, col);
00099       return coeffRef(row,col);
00100     }
00101 
00102     #ifndef EIGEN_PARSED_BY_DOXYGEN
00103     EIGEN_DEVICE_FUNC
00104     inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
00105     EIGEN_DEVICE_FUNC
00106     inline Derived& derived() { return *static_cast<Derived*>(this); }
00107     #endif // not EIGEN_PARSED_BY_DOXYGEN
00108 
00109     template<typename DenseDerived>
00110     EIGEN_DEVICE_FUNC
00111     void evalTo(MatrixBase<DenseDerived> &other) const;
00112     template<typename DenseDerived>
00113     EIGEN_DEVICE_FUNC
00114     void evalToLazy(MatrixBase<DenseDerived> &other) const;
00115 
00116     EIGEN_DEVICE_FUNC
00117     DenseMatrixType toDenseMatrix() const
00118     {
00119       DenseMatrixType res(rows(), cols());
00120       evalToLazy(res);
00121       return res;
00122     }
00123 
00124   protected:
00125 
00126     void check_coordinates(Index row, Index col) const
00127     {
00128       EIGEN_ONLY_USED_FOR_DEBUG(row);
00129       EIGEN_ONLY_USED_FOR_DEBUG(col);
00130       eigen_assert(col>=0 && col<cols() && row>=0 && row<rows());
00131       const int mode = int(Mode) & ~SelfAdjoint;
00132       EIGEN_ONLY_USED_FOR_DEBUG(mode);
00133       eigen_assert((mode==Upper && col>=row)
00134                 || (mode==Lower && col<=row)
00135                 || ((mode==StrictlyUpper || mode==UnitUpper) && col>row)
00136                 || ((mode==StrictlyLower || mode==UnitLower) && col<row));
00137     }
00138 
00139     #ifdef EIGEN_INTERNAL_DEBUGGING
00140     void check_coordinates_internal(Index row, Index col) const
00141     {
00142       check_coordinates(row, col);
00143     }
00144     #else
00145     void check_coordinates_internal(Index , Index ) const {}
00146     #endif
00147 
00148 };
00149 
00167 namespace internal {
00168 template<typename MatrixType, unsigned int _Mode>
00169 struct traits<TriangularView<MatrixType, _Mode> > : traits<MatrixType>
00170 {
00171   typedef typename ref_selector<MatrixType>::non_const_type MatrixTypeNested;
00172   typedef typename remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef;
00173   typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
00174   typedef typename MatrixType::PlainObject FullMatrixType;
00175   typedef MatrixType ExpressionType;
00176   enum {
00177     Mode = _Mode,
00178     FlagsLvalueBit = is_lvalue<MatrixType>::value ? LvalueBit : 0,
00179     Flags = (MatrixTypeNestedCleaned::Flags & (HereditaryBits | FlagsLvalueBit) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)))
00180   };
00181 };
00182 }
00183 
00184 template<typename _MatrixType, unsigned int _Mode, typename StorageKind> class TriangularViewImpl;
00185 
00186 template<typename _MatrixType, unsigned int _Mode> class TriangularView
00187   : public TriangularViewImpl<_MatrixType, _Mode, typename internal::traits<_MatrixType>::StorageKind >
00188 {
00189   public:
00190 
00191     typedef TriangularViewImpl<_MatrixType, _Mode, typename internal::traits<_MatrixType>::StorageKind > Base;
00192     typedef typename internal::traits<TriangularView>::Scalar Scalar;
00193     typedef _MatrixType MatrixType;
00194 
00195   protected:
00196     typedef typename internal::traits<TriangularView>::MatrixTypeNested MatrixTypeNested;
00197     typedef typename internal::traits<TriangularView>::MatrixTypeNestedNonRef MatrixTypeNestedNonRef;
00198 
00199     typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType;
00200     
00201   public:
00202 
00203     typedef typename internal::traits<TriangularView>::StorageKind StorageKind;
00204     typedef typename internal::traits<TriangularView>::MatrixTypeNestedCleaned NestedExpression;
00205 
00206     enum {
00207       Mode = _Mode,
00208       Flags = internal::traits<TriangularView>::Flags,
00209       TransposeMode = (Mode & Upper ? Lower : 0)
00210                     | (Mode & Lower ? Upper : 0)
00211                     | (Mode & (UnitDiag))
00212                     | (Mode & (ZeroDiag)),
00213       IsVectorAtCompileTime = false
00214     };
00215 
00216     EIGEN_DEVICE_FUNC
00217     explicit inline TriangularView(MatrixType& matrix) : m_matrix(matrix)
00218     {}
00219     
00220     using Base::operator=;
00221     TriangularView& operator=(const TriangularView &other)
00222     { return Base::operator=(other); }
00223 
00225     EIGEN_DEVICE_FUNC
00226     inline Index rows() const { return m_matrix.rows(); }
00228     EIGEN_DEVICE_FUNC
00229     inline Index cols() const { return m_matrix.cols(); }
00230 
00232     EIGEN_DEVICE_FUNC
00233     const NestedExpression& nestedExpression() const { return m_matrix; }
00234 
00236     EIGEN_DEVICE_FUNC
00237     NestedExpression& nestedExpression() { return m_matrix; }
00238     
00239     typedef TriangularView<const MatrixConjugateReturnType,Mode> ConjugateReturnType;
00241     EIGEN_DEVICE_FUNC
00242     inline const ConjugateReturnType conjugate() const
00243     { return ConjugateReturnType(m_matrix.conjugate()); }
00244 
00245     typedef TriangularView<const typename MatrixType::AdjointReturnType,TransposeMode> AdjointReturnType;
00247     EIGEN_DEVICE_FUNC
00248     inline const AdjointReturnType adjoint() const
00249     { return AdjointReturnType(m_matrix.adjoint()); }
00250 
00251     typedef TriangularView<typename MatrixType::TransposeReturnType,TransposeMode> TransposeReturnType;
00253     EIGEN_DEVICE_FUNC
00254     inline TransposeReturnType transpose()
00255     {
00256       EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
00257       typename MatrixType::TransposeReturnType tmp(m_matrix);
00258       return TransposeReturnType(tmp);
00259     }
00260     
00261     typedef TriangularView<const typename MatrixType::ConstTransposeReturnType,TransposeMode> ConstTransposeReturnType;
00263     EIGEN_DEVICE_FUNC
00264     inline const ConstTransposeReturnType transpose() const
00265     {
00266       return ConstTransposeReturnType(m_matrix.transpose());
00267     }
00268 
00269     template<typename Other>
00270     EIGEN_DEVICE_FUNC
00271     inline const Solve<TriangularView, Other> 
00272     solve(const MatrixBase<Other>& other) const
00273     { return Solve<TriangularView, Other>(*this, other.derived()); }
00274     
00275   // workaround MSVC ICE
00276   #if EIGEN_COMP_MSVC
00277     template<int Side, typename Other>
00278     EIGEN_DEVICE_FUNC
00279     inline const internal::triangular_solve_retval<Side,TriangularView, Other>
00280     solve(const MatrixBase<Other>& other) const
00281     { return Base::template solve<Side>(other); }
00282   #else
00283     using Base::solve;
00284   #endif
00285 
00290     EIGEN_DEVICE_FUNC
00291     SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView()
00292     {
00293       EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
00294       return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
00295     }
00296 
00298     EIGEN_DEVICE_FUNC
00299     const SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() const
00300     {
00301       EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
00302       return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
00303     }
00304 
00305 
00308     EIGEN_DEVICE_FUNC
00309     Scalar determinant() const
00310     {
00311       if (Mode & UnitDiag)
00312         return 1;
00313       else if (Mode & ZeroDiag)
00314         return 0;
00315       else
00316         return m_matrix.diagonal().prod();
00317     }
00318       
00319   protected:
00320 
00321     MatrixTypeNested m_matrix;
00322 };
00323 
00333 template<typename _MatrixType, unsigned int _Mode> class TriangularViewImpl<_MatrixType,_Mode,Dense>
00334   : public TriangularBase<TriangularView<_MatrixType, _Mode> >
00335 {
00336   public:
00337 
00338     typedef TriangularView<_MatrixType, _Mode> TriangularViewType;
00339     typedef TriangularBase<TriangularViewType> Base;
00340     typedef typename internal::traits<TriangularViewType>::Scalar Scalar;
00341 
00342     typedef _MatrixType MatrixType;
00343     typedef typename MatrixType::PlainObject DenseMatrixType;
00344     typedef DenseMatrixType PlainObject;
00345 
00346   public:
00347     using Base::evalToLazy;
00348     using Base::derived;
00349 
00350     typedef typename internal::traits<TriangularViewType>::StorageKind StorageKind;
00351 
00352     enum {
00353       Mode = _Mode,
00354       Flags = internal::traits<TriangularViewType>::Flags
00355     };
00356 
00359     EIGEN_DEVICE_FUNC
00360     inline Index outerStride() const { return derived().nestedExpression().outerStride(); }
00363     EIGEN_DEVICE_FUNC
00364     inline Index innerStride() const { return derived().nestedExpression().innerStride(); }
00365 
00367     template<typename Other>
00368     EIGEN_DEVICE_FUNC
00369     TriangularViewType&  operator+=(const DenseBase<Other>& other) {
00370       internal::call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op<Scalar,typename Other::Scalar>());
00371       return derived();
00372     }
00374     template<typename Other>
00375     EIGEN_DEVICE_FUNC
00376     TriangularViewType&  operator-=(const DenseBase<Other>& other) {
00377       internal::call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op<Scalar,typename Other::Scalar>());
00378       return derived();
00379     }
00380     
00382     EIGEN_DEVICE_FUNC
00383     TriangularViewType&  operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() * other; }
00385     EIGEN_DEVICE_FUNC
00386     TriangularViewType&  operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() / other; }
00387 
00389     EIGEN_DEVICE_FUNC
00390     void fill(const Scalar& value) { setConstant(value); }
00392     EIGEN_DEVICE_FUNC
00393     TriangularViewType& setConstant(const Scalar& value)
00394     { return *this = MatrixType::Constant(derived().rows(), derived().cols(), value); }
00396     EIGEN_DEVICE_FUNC
00397     TriangularViewType& setZero() { return setConstant(Scalar(0)); }
00399     EIGEN_DEVICE_FUNC
00400     TriangularViewType& setOnes() { return setConstant(Scalar(1)); }
00401 
00405     EIGEN_DEVICE_FUNC
00406     inline Scalar coeff(Index row, Index col) const
00407     {
00408       Base::check_coordinates_internal(row, col);
00409       return derived().nestedExpression().coeff(row, col);
00410     }
00411 
00415     EIGEN_DEVICE_FUNC
00416     inline Scalar& coeffRef(Index row, Index col)
00417     {
00418       EIGEN_STATIC_ASSERT_LVALUE(TriangularViewType);
00419       Base::check_coordinates_internal(row, col);
00420       return derived().nestedExpression().coeffRef(row, col);
00421     }
00422 
00424     template<typename OtherDerived>
00425     EIGEN_DEVICE_FUNC
00426     TriangularViewType& operator=(const TriangularBase<OtherDerived>& other);
00427 
00429     template<typename OtherDerived>
00430     EIGEN_DEVICE_FUNC
00431     TriangularViewType& operator=(const MatrixBase<OtherDerived>& other);
00432 
00433 #ifndef EIGEN_PARSED_BY_DOXYGEN
00434     EIGEN_DEVICE_FUNC
00435     TriangularViewType& operator=(const TriangularViewImpl& other)
00436     { return *this = other.derived().nestedExpression(); }
00437 
00439     template<typename OtherDerived>
00440     EIGEN_DEVICE_FUNC
00441     void lazyAssign(const TriangularBase<OtherDerived>& other);
00442 
00444     template<typename OtherDerived>
00445     EIGEN_DEVICE_FUNC
00446     void lazyAssign(const MatrixBase<OtherDerived>& other);
00447 #endif
00448 
00450     template<typename OtherDerived>
00451     EIGEN_DEVICE_FUNC
00452     const Product<TriangularViewType,OtherDerived>
00453     operator*(const MatrixBase<OtherDerived>& rhs) const
00454     {
00455       return Product<TriangularViewType,OtherDerived>(derived(), rhs.derived());
00456     }
00457 
00459     template<typename OtherDerived> friend
00460     EIGEN_DEVICE_FUNC
00461     const Product<OtherDerived,TriangularViewType>
00462     operator*(const MatrixBase<OtherDerived>& lhs, const TriangularViewImpl& rhs)
00463     {
00464       return Product<OtherDerived,TriangularViewType>(lhs.derived(),rhs.derived());
00465     }
00466 
00490     template<int Side, typename Other>
00491     EIGEN_DEVICE_FUNC
00492     inline const internal::triangular_solve_retval<Side,TriangularViewType, Other>
00493     solve(const MatrixBase<Other>& other) const;
00494 
00504     template<int Side, typename OtherDerived>
00505     EIGEN_DEVICE_FUNC
00506     void solveInPlace(const MatrixBase<OtherDerived>& other) const;
00507 
00508     template<typename OtherDerived>
00509     EIGEN_DEVICE_FUNC
00510     void solveInPlace(const MatrixBase<OtherDerived>& other) const
00511     { return solveInPlace<OnTheLeft>(other); }
00512 
00514     template<typename OtherDerived>
00515     EIGEN_DEVICE_FUNC
00516 #ifdef EIGEN_PARSED_BY_DOXYGEN
00517     void swap(TriangularBase<OtherDerived> &other)
00518 #else
00519     void swap(TriangularBase<OtherDerived> const & other)
00520 #endif
00521     {
00522       EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
00523       call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
00524     }
00525 
00528     template<typename OtherDerived>
00529     EIGEN_DEVICE_FUNC
00530     void swap(MatrixBase<OtherDerived> const & other)
00531     {
00532       EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
00533       call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
00534     }
00535 
00536     template<typename RhsType, typename DstType>
00537     EIGEN_DEVICE_FUNC
00538     EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const {
00539       if(!internal::is_same_dense(dst,rhs))
00540         dst = rhs;
00541       this->solveInPlace(dst);
00542     }
00543 
00544     template<typename ProductType>
00545     EIGEN_DEVICE_FUNC
00546     EIGEN_STRONG_INLINE TriangularViewType& _assignProduct(const ProductType& prod, const Scalar& alpha, bool beta);
00547 };
00548 
00549 /***************************************************************************
00550 * Implementation of triangular evaluation/assignment
00551 ***************************************************************************/
00552 
00553 #ifndef EIGEN_PARSED_BY_DOXYGEN
00554 // FIXME should we keep that possibility
00555 template<typename MatrixType, unsigned int Mode>
00556 template<typename OtherDerived>
00557 inline TriangularView<MatrixType, Mode>&
00558 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const MatrixBase<OtherDerived>& other)
00559 {
00560   internal::call_assignment_no_alias(derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
00561   return derived();
00562 }
00563 
00564 // FIXME should we keep that possibility
00565 template<typename MatrixType, unsigned int Mode>
00566 template<typename OtherDerived>
00567 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const MatrixBase<OtherDerived>& other)
00568 {
00569   internal::call_assignment_no_alias(derived(), other.template triangularView<Mode>());
00570 }
00571 
00572 
00573 
00574 template<typename MatrixType, unsigned int Mode>
00575 template<typename OtherDerived>
00576 inline TriangularView<MatrixType, Mode>&
00577 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const TriangularBase<OtherDerived>& other)
00578 {
00579   eigen_assert(Mode == int(OtherDerived::Mode));
00580   internal::call_assignment(derived(), other.derived());
00581   return derived();
00582 }
00583 
00584 template<typename MatrixType, unsigned int Mode>
00585 template<typename OtherDerived>
00586 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const TriangularBase<OtherDerived>& other)
00587 {
00588   eigen_assert(Mode == int(OtherDerived::Mode));
00589   internal::call_assignment_no_alias(derived(), other.derived());
00590 }
00591 #endif
00592 
00593 /***************************************************************************
00594 * Implementation of TriangularBase methods
00595 ***************************************************************************/
00596 
00599 template<typename Derived>
00600 template<typename DenseDerived>
00601 void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
00602 {
00603   evalToLazy(other.derived());
00604 }
00605 
00606 /***************************************************************************
00607 * Implementation of TriangularView methods
00608 ***************************************************************************/
00609 
00610 /***************************************************************************
00611 * Implementation of MatrixBase methods
00612 ***************************************************************************/
00613 
00625 template<typename Derived>
00626 template<unsigned int Mode>
00627 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
00628 MatrixBase<Derived>::triangularView()
00629 {
00630   return typename TriangularViewReturnType<Mode>::Type(derived());
00631 }
00632 
00634 template<typename Derived>
00635 template<unsigned int Mode>
00636 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
00637 MatrixBase<Derived>::triangularView() const
00638 {
00639   return typename ConstTriangularViewReturnType<Mode>::Type(derived());
00640 }
00641 
00647 template<typename Derived>
00648 bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const
00649 {
00650   RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
00651   for(Index j = 0; j < cols(); ++j)
00652   {
00653     Index maxi = numext::mini(j, rows()-1);
00654     for(Index i = 0; i <= maxi; ++i)
00655     {
00656       RealScalar absValue = numext::abs(coeff(i,j));
00657       if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
00658     }
00659   }
00660   RealScalar threshold = maxAbsOnUpperPart * prec;
00661   for(Index j = 0; j < cols(); ++j)
00662     for(Index i = j+1; i < rows(); ++i)
00663       if(numext::abs(coeff(i, j)) > threshold) return false;
00664   return true;
00665 }
00666 
00672 template<typename Derived>
00673 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const
00674 {
00675   RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
00676   for(Index j = 0; j < cols(); ++j)
00677     for(Index i = j; i < rows(); ++i)
00678     {
00679       RealScalar absValue = numext::abs(coeff(i,j));
00680       if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
00681     }
00682   RealScalar threshold = maxAbsOnLowerPart * prec;
00683   for(Index j = 1; j < cols(); ++j)
00684   {
00685     Index maxi = numext::mini(j, rows()-1);
00686     for(Index i = 0; i < maxi; ++i)
00687       if(numext::abs(coeff(i, j)) > threshold) return false;
00688   }
00689   return true;
00690 }
00691 
00692 
00693 /***************************************************************************
00694 ****************************************************************************
00695 * Evaluators and Assignment of triangular expressions
00696 ***************************************************************************
00697 ***************************************************************************/
00698 
00699 namespace internal {
00700 
00701   
00702 // TODO currently a triangular expression has the form TriangularView<.,.>
00703 //      in the future triangular-ness should be defined by the expression traits
00704 //      such that Transpose<TriangularView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
00705 template<typename MatrixType, unsigned int Mode>
00706 struct evaluator_traits<TriangularView<MatrixType,Mode> >
00707 {
00708   typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
00709   typedef typename glue_shapes<typename evaluator_traits<MatrixType>::Shape, TriangularShape>::type Shape;
00710 };
00711 
00712 template<typename MatrixType, unsigned int Mode>
00713 struct unary_evaluator<TriangularView<MatrixType,Mode>, IndexBased>
00714  : evaluator<typename internal::remove_all<MatrixType>::type>
00715 {
00716   typedef TriangularView<MatrixType,Mode> XprType;
00717   typedef evaluator<typename internal::remove_all<MatrixType>::type> Base;
00718   unary_evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {}
00719 };
00720 
00721 // Additional assignment kinds:
00722 struct Triangular2Triangular    {};
00723 struct Triangular2Dense         {};
00724 struct Dense2Triangular         {};
00725 
00726 
00727 template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop;
00728 
00729  
00735 template<int UpLo, int Mode, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version = Specialized>
00736 class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version>
00737 {
00738 protected:
00739   typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base;
00740   typedef typename Base::DstXprType DstXprType;
00741   typedef typename Base::SrcXprType SrcXprType;
00742   using Base::m_dst;
00743   using Base::m_src;
00744   using Base::m_functor;
00745 public:
00746   
00747   typedef typename Base::DstEvaluatorType DstEvaluatorType;
00748   typedef typename Base::SrcEvaluatorType SrcEvaluatorType;
00749   typedef typename Base::Scalar Scalar;
00750   typedef typename Base::AssignmentTraits AssignmentTraits;
00751   
00752   
00753   EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr)
00754     : Base(dst, src, func, dstExpr)
00755   {}
00756   
00757 #ifdef EIGEN_INTERNAL_DEBUGGING
00758   EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col)
00759   {
00760     eigen_internal_assert(row!=col);
00761     Base::assignCoeff(row,col);
00762   }
00763 #else
00764   using Base::assignCoeff;
00765 #endif
00766   
00767   EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id)
00768   {
00769          if(Mode==UnitDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(1));
00770     else if(Mode==ZeroDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(0));
00771     else if(Mode==0)                       Base::assignCoeff(id,id);
00772   }
00773   
00774   EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index row, Index col)
00775   { 
00776     eigen_internal_assert(row!=col);
00777     if(SetOpposite)
00778       m_functor.assignCoeff(m_dst.coeffRef(row,col), Scalar(0));
00779   }
00780 };
00781 
00782 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType, typename Functor>
00783 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
00784 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src, const Functor &func)
00785 {
00786   typedef evaluator<DstXprType> DstEvaluatorType;
00787   typedef evaluator<SrcXprType> SrcEvaluatorType;
00788 
00789   SrcEvaluatorType srcEvaluator(src);
00790 
00791   Index dstRows = src.rows();
00792   Index dstCols = src.cols();
00793   if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
00794     dst.resize(dstRows, dstCols);
00795   DstEvaluatorType dstEvaluator(dst);
00796     
00797   typedef triangular_dense_assignment_kernel< Mode&(Lower|Upper),Mode&(UnitDiag|ZeroDiag|SelfAdjoint),SetOpposite,
00798                                               DstEvaluatorType,SrcEvaluatorType,Functor> Kernel;
00799   Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived());
00800   
00801   enum {
00802       unroll = DstXprType::SizeAtCompileTime != Dynamic
00803             && SrcEvaluatorType::CoeffReadCost < HugeCost
00804             && DstXprType::SizeAtCompileTime * (DstEvaluatorType::CoeffReadCost+SrcEvaluatorType::CoeffReadCost) / 2 <= EIGEN_UNROLLING_LIMIT
00805     };
00806   
00807   triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, SetOpposite>::run(kernel);
00808 }
00809 
00810 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType>
00811 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
00812 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src)
00813 {
00814   call_triangular_assignment_loop<Mode,SetOpposite>(dst, src, internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>());
00815 }
00816 
00817 template<> struct AssignmentKind<TriangularShape,TriangularShape> { typedef Triangular2Triangular Kind; };
00818 template<> struct AssignmentKind<DenseShape,TriangularShape>      { typedef Triangular2Dense      Kind; };
00819 template<> struct AssignmentKind<TriangularShape,DenseShape>      { typedef Dense2Triangular      Kind; };
00820 
00821 
00822 template< typename DstXprType, typename SrcXprType, typename Functor>
00823 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Triangular>
00824 {
00825   EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
00826   {
00827     eigen_assert(int(DstXprType::Mode) == int(SrcXprType::Mode));
00828     
00829     call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);  
00830   }
00831 };
00832 
00833 template< typename DstXprType, typename SrcXprType, typename Functor>
00834 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Dense>
00835 {
00836   EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
00837   {
00838     call_triangular_assignment_loop<SrcXprType::Mode, (SrcXprType::Mode&SelfAdjoint)==0>(dst, src, func);  
00839   }
00840 };
00841 
00842 template< typename DstXprType, typename SrcXprType, typename Functor>
00843 struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular>
00844 {
00845   EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
00846   {
00847     call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);  
00848   }
00849 };
00850 
00851 
00852 template<typename Kernel, unsigned int Mode, int UnrollCount, bool SetOpposite>
00853 struct triangular_assignment_loop
00854 {
00855   // FIXME: this is not very clean, perhaps this information should be provided by the kernel?
00856   typedef typename Kernel::DstEvaluatorType DstEvaluatorType;
00857   typedef typename DstEvaluatorType::XprType DstXprType;
00858   
00859   enum {
00860     col = (UnrollCount-1) / DstXprType::RowsAtCompileTime,
00861     row = (UnrollCount-1) % DstXprType::RowsAtCompileTime
00862   };
00863   
00864   typedef typename Kernel::Scalar Scalar;
00865 
00866   EIGEN_DEVICE_FUNC
00867   static inline void run(Kernel &kernel)
00868   {
00869     triangular_assignment_loop<Kernel, Mode, UnrollCount-1, SetOpposite>::run(kernel);
00870     
00871     if(row==col)
00872       kernel.assignDiagonalCoeff(row);
00873     else if( ((Mode&Lower) && row>col) || ((Mode&Upper) && row<col) )
00874       kernel.assignCoeff(row,col);
00875     else if(SetOpposite)
00876       kernel.assignOppositeCoeff(row,col);
00877   }
00878 };
00879 
00880 // prevent buggy user code from causing an infinite recursion
00881 template<typename Kernel, unsigned int Mode, bool SetOpposite>
00882 struct triangular_assignment_loop<Kernel, Mode, 0, SetOpposite>
00883 {
00884   EIGEN_DEVICE_FUNC
00885   static inline void run(Kernel &) {}
00886 };
00887 
00888 
00889 
00890 // TODO: experiment with a recursive assignment procedure splitting the current
00891 //       triangular part into one rectangular and two triangular parts.
00892 
00893 
00894 template<typename Kernel, unsigned int Mode, bool SetOpposite>
00895 struct triangular_assignment_loop<Kernel, Mode, Dynamic, SetOpposite>
00896 {
00897   typedef typename Kernel::Scalar Scalar;
00898   EIGEN_DEVICE_FUNC
00899   static inline void run(Kernel &kernel)
00900   {
00901     for(Index j = 0; j < kernel.cols(); ++j)
00902     {
00903       Index maxi = numext::mini(j, kernel.rows());
00904       Index i = 0;
00905       if (((Mode&Lower) && SetOpposite) || (Mode&Upper))
00906       {
00907         for(; i < maxi; ++i)
00908           if(Mode&Upper) kernel.assignCoeff(i, j);
00909           else           kernel.assignOppositeCoeff(i, j);
00910       }
00911       else
00912         i = maxi;
00913       
00914       if(i<kernel.rows()) // then i==j
00915         kernel.assignDiagonalCoeff(i++);
00916       
00917       if (((Mode&Upper) && SetOpposite) || (Mode&Lower))
00918       {
00919         for(; i < kernel.rows(); ++i)
00920           if(Mode&Lower) kernel.assignCoeff(i, j);
00921           else           kernel.assignOppositeCoeff(i, j);
00922       }
00923     }
00924   }
00925 };
00926 
00927 } // end namespace internal
00928 
00931 template<typename Derived>
00932 template<typename DenseDerived>
00933 void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
00934 {
00935   other.derived().resize(this->rows(), this->cols());
00936   internal::call_triangular_assignment_loop<Derived::Mode,(Derived::Mode&SelfAdjoint)==0 /* SetOpposite */>(other.derived(), derived().nestedExpression());
00937 }
00938 
00939 namespace internal {
00940   
00941 // Triangular = Product
00942 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
00943 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
00944 {
00945   typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
00946   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,typename SrcXprType::Scalar> &)
00947   {
00948     Index dstRows = src.rows();
00949     Index dstCols = src.cols();
00950     if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
00951       dst.resize(dstRows, dstCols);
00952 
00953     dst._assignProduct(src, 1, 0);
00954   }
00955 };
00956 
00957 // Triangular += Product
00958 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
00959 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::add_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
00960 {
00961   typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
00962   static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar,typename SrcXprType::Scalar> &)
00963   {
00964     dst._assignProduct(src, 1, 1);
00965   }
00966 };
00967 
00968 // Triangular -= Product
00969 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
00970 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::sub_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
00971 {
00972   typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
00973   static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,typename SrcXprType::Scalar> &)
00974   {
00975     dst._assignProduct(src, -1, 1);
00976   }
00977 };
00978 
00979 } // end namespace internal
00980 
00981 } // end namespace Eigen
00982 
00983 #endif // EIGEN_TRIANGULARMATRIX_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends