![]() |
Eigen
3.3.3
|
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