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