1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008 Benoit Jacob <jacob.benoit.1 (at) gmail.com> 5 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud (at) inria.fr> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_TRIANGULARMATRIX_H 12 #define EIGEN_TRIANGULARMATRIX_H 13 14 namespace Eigen { 15 16 namespace internal { 17 18 template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval; 19 20 } 21 22 /** \internal 23 * 24 * \class TriangularBase 25 * \ingroup Core_Module 26 * 27 * \brief Base class for triangular part in a matrix 28 */ 29 template<typename Derived> class TriangularBase : public EigenBase<Derived> 30 { 31 public: 32 33 enum { 34 Mode = internal::traits<Derived>::Mode, 35 CoeffReadCost = internal::traits<Derived>::CoeffReadCost, 36 RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime, 37 ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime, 38 MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime, 39 MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime 40 }; 41 typedef typename internal::traits<Derived>::Scalar Scalar; 42 typedef typename internal::traits<Derived>::StorageKind StorageKind; 43 typedef typename internal::traits<Derived>::Index Index; 44 typedef typename internal::traits<Derived>::DenseMatrixType DenseMatrixType; 45 typedef DenseMatrixType DenseType; 46 47 inline TriangularBase() { eigen_assert(!((Mode&UnitDiag) && (Mode&ZeroDiag))); } 48 49 inline Index rows() const { return derived().rows(); } 50 inline Index cols() const { return derived().cols(); } 51 inline Index outerStride() const { return derived().outerStride(); } 52 inline Index innerStride() const { return derived().innerStride(); } 53 54 inline Scalar coeff(Index row, Index col) const { return derived().coeff(row,col); } 55 inline Scalar& coeffRef(Index row, Index col) { return derived().coeffRef(row,col); } 56 57 /** \see MatrixBase::copyCoeff(row,col) 58 */ 59 template<typename Other> 60 EIGEN_STRONG_INLINE void copyCoeff(Index row, Index col, Other& other) 61 { 62 derived().coeffRef(row, col) = other.coeff(row, col); 63 } 64 65 inline Scalar operator()(Index row, Index col) const 66 { 67 check_coordinates(row, col); 68 return coeff(row,col); 69 } 70 inline Scalar& operator()(Index row, Index col) 71 { 72 check_coordinates(row, col); 73 return coeffRef(row,col); 74 } 75 76 #ifndef EIGEN_PARSED_BY_DOXYGEN 77 inline const Derived& derived() const { return *static_cast<const Derived*>(this); } 78 inline Derived& derived() { return *static_cast<Derived*>(this); } 79 #endif // not EIGEN_PARSED_BY_DOXYGEN 80 81 template<typename DenseDerived> 82 void evalTo(MatrixBase<DenseDerived> &other) const; 83 template<typename DenseDerived> 84 void evalToLazy(MatrixBase<DenseDerived> &other) const; 85 86 DenseMatrixType toDenseMatrix() const 87 { 88 DenseMatrixType res(rows(), cols()); 89 evalToLazy(res); 90 return res; 91 } 92 93 protected: 94 95 void check_coordinates(Index row, Index col) const 96 { 97 EIGEN_ONLY_USED_FOR_DEBUG(row); 98 EIGEN_ONLY_USED_FOR_DEBUG(col); 99 eigen_assert(col>=0 && col<cols() && row>=0 && row<rows()); 100 const int mode = int(Mode) & ~SelfAdjoint; 101 EIGEN_ONLY_USED_FOR_DEBUG(mode); 102 eigen_assert((mode==Upper && col>=row) 103 || (mode==Lower && col<=row) 104 || ((mode==StrictlyUpper || mode==UnitUpper) && col>row) 105 || ((mode==StrictlyLower || mode==UnitLower) && col<row)); 106 } 107 108 #ifdef EIGEN_INTERNAL_DEBUGGING 109 void check_coordinates_internal(Index row, Index col) const 110 { 111 check_coordinates(row, col); 112 } 113 #else 114 void check_coordinates_internal(Index , Index ) const {} 115 #endif 116 117 }; 118 119 /** \class TriangularView 120 * \ingroup Core_Module 121 * 122 * \brief Base class for triangular part in a matrix 123 * 124 * \param MatrixType the type of the object in which we are taking the triangular part 125 * \param Mode the kind of triangular matrix expression to construct. Can be #Upper, 126 * #Lower, #UnitUpper, #UnitLower, #StrictlyUpper, or #StrictlyLower. 127 * This is in fact a bit field; it must have either #Upper or #Lower, 128 * and additionnaly it may have #UnitDiag or #ZeroDiag or neither. 129 * 130 * This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular 131 * matrices one should speak of "trapezoid" parts. This class is the return type 132 * of MatrixBase::triangularView() and most of the time this is the only way it is used. 133 * 134 * \sa MatrixBase::triangularView() 135 */ 136 namespace internal { 137 template<typename MatrixType, unsigned int _Mode> 138 struct traits<TriangularView<MatrixType, _Mode> > : traits<MatrixType> 139 { 140 typedef typename nested<MatrixType>::type MatrixTypeNested; 141 typedef typename remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef; 142 typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned; 143 typedef MatrixType ExpressionType; 144 typedef typename MatrixType::PlainObject DenseMatrixType; 145 enum { 146 Mode = _Mode, 147 Flags = (MatrixTypeNestedCleaned::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) | Mode, 148 CoeffReadCost = MatrixTypeNestedCleaned::CoeffReadCost 149 }; 150 }; 151 } 152 153 template<int Mode, bool LhsIsTriangular, 154 typename Lhs, bool LhsIsVector, 155 typename Rhs, bool RhsIsVector> 156 struct TriangularProduct; 157 158 template<typename _MatrixType, unsigned int _Mode> class TriangularView 159 : public TriangularBase<TriangularView<_MatrixType, _Mode> > 160 { 161 public: 162 163 typedef TriangularBase<TriangularView> Base; 164 typedef typename internal::traits<TriangularView>::Scalar Scalar; 165 166 typedef _MatrixType MatrixType; 167 typedef typename internal::traits<TriangularView>::DenseMatrixType DenseMatrixType; 168 typedef DenseMatrixType PlainObject; 169 170 protected: 171 typedef typename internal::traits<TriangularView>::MatrixTypeNested MatrixTypeNested; 172 typedef typename internal::traits<TriangularView>::MatrixTypeNestedNonRef MatrixTypeNestedNonRef; 173 typedef typename internal::traits<TriangularView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned; 174 175 typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType; 176 177 public: 178 using Base::evalToLazy; 179 180 181 typedef typename internal::traits<TriangularView>::StorageKind StorageKind; 182 typedef typename internal::traits<TriangularView>::Index Index; 183 184 enum { 185 Mode = _Mode, 186 TransposeMode = (Mode & Upper ? Lower : 0) 187 | (Mode & Lower ? Upper : 0) 188 | (Mode & (UnitDiag)) 189 | (Mode & (ZeroDiag)) 190 }; 191 192 inline TriangularView(const MatrixType& matrix) : m_matrix(matrix) 193 {} 194 195 inline Index rows() const { return m_matrix.rows(); } 196 inline Index cols() const { return m_matrix.cols(); } 197 inline Index outerStride() const { return m_matrix.outerStride(); } 198 inline Index innerStride() const { return m_matrix.innerStride(); } 199 200 /** \sa MatrixBase::operator+=() */ 201 template<typename Other> TriangularView& operator+=(const DenseBase<Other>& other) { return *this = m_matrix + other.derived(); } 202 /** \sa MatrixBase::operator-=() */ 203 template<typename Other> TriangularView& operator-=(const DenseBase<Other>& other) { return *this = m_matrix - other.derived(); } 204 /** \sa MatrixBase::operator*=() */ 205 TriangularView& operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix * other; } 206 /** \sa MatrixBase::operator/=() */ 207 TriangularView& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix / other; } 208 209 /** \sa MatrixBase::fill() */ 210 void fill(const Scalar& value) { setConstant(value); } 211 /** \sa MatrixBase::setConstant() */ 212 TriangularView& setConstant(const Scalar& value) 213 { return *this = MatrixType::Constant(rows(), cols(), value); } 214 /** \sa MatrixBase::setZero() */ 215 TriangularView& setZero() { return setConstant(Scalar(0)); } 216 /** \sa MatrixBase::setOnes() */ 217 TriangularView& setOnes() { return setConstant(Scalar(1)); } 218 219 /** \sa MatrixBase::coeff() 220 * \warning the coordinates must fit into the referenced triangular part 221 */ 222 inline Scalar coeff(Index row, Index col) const 223 { 224 Base::check_coordinates_internal(row, col); 225 return m_matrix.coeff(row, col); 226 } 227 228 /** \sa MatrixBase::coeffRef() 229 * \warning the coordinates must fit into the referenced triangular part 230 */ 231 inline Scalar& coeffRef(Index row, Index col) 232 { 233 Base::check_coordinates_internal(row, col); 234 return m_matrix.const_cast_derived().coeffRef(row, col); 235 } 236 237 const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; } 238 MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); } 239 240 /** Assigns a triangular matrix to a triangular part of a dense matrix */ 241 template<typename OtherDerived> 242 TriangularView& operator=(const TriangularBase<OtherDerived>& other); 243 244 template<typename OtherDerived> 245 TriangularView& operator=(const MatrixBase<OtherDerived>& other); 246 247 TriangularView& operator=(const TriangularView& other) 248 { return *this = other.nestedExpression(); } 249 250 template<typename OtherDerived> 251 void lazyAssign(const TriangularBase<OtherDerived>& other); 252 253 template<typename OtherDerived> 254 void lazyAssign(const MatrixBase<OtherDerived>& other); 255 256 /** \sa MatrixBase::conjugate() */ 257 inline TriangularView<MatrixConjugateReturnType,Mode> conjugate() 258 { return m_matrix.conjugate(); } 259 /** \sa MatrixBase::conjugate() const */ 260 inline const TriangularView<MatrixConjugateReturnType,Mode> conjugate() const 261 { return m_matrix.conjugate(); } 262 263 /** \sa MatrixBase::adjoint() const */ 264 inline const TriangularView<const typename MatrixType::AdjointReturnType,TransposeMode> adjoint() const 265 { return m_matrix.adjoint(); } 266 267 /** \sa MatrixBase::transpose() */ 268 inline TriangularView<Transpose<MatrixType>,TransposeMode> transpose() 269 { 270 EIGEN_STATIC_ASSERT_LVALUE(MatrixType) 271 return m_matrix.const_cast_derived().transpose(); 272 } 273 /** \sa MatrixBase::transpose() const */ 274 inline const TriangularView<Transpose<MatrixType>,TransposeMode> transpose() const 275 { 276 return m_matrix.transpose(); 277 } 278 279 /** Efficient triangular matrix times vector/matrix product */ 280 template<typename OtherDerived> 281 TriangularProduct<Mode, true, MatrixType, false, OtherDerived, OtherDerived::ColsAtCompileTime==1> 282 operator*(const MatrixBase<OtherDerived>& rhs) const 283 { 284 return TriangularProduct 285 <Mode, true, MatrixType, false, OtherDerived, OtherDerived::ColsAtCompileTime==1> 286 (m_matrix, rhs.derived()); 287 } 288 289 /** Efficient vector/matrix times triangular matrix product */ 290 template<typename OtherDerived> friend 291 TriangularProduct<Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false> 292 operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs) 293 { 294 return TriangularProduct 295 <Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false> 296 (lhs.derived(),rhs.m_matrix); 297 } 298 299 #ifdef EIGEN2_SUPPORT 300 template<typename OtherDerived> 301 struct eigen2_product_return_type 302 { 303 typedef typename TriangularView<MatrixType,Mode>::DenseMatrixType DenseMatrixType; 304 typedef typename OtherDerived::PlainObject::DenseType OtherPlainObject; 305 typedef typename ProductReturnType<DenseMatrixType, OtherPlainObject>::Type ProdRetType; 306 typedef typename ProdRetType::PlainObject type; 307 }; 308 template<typename OtherDerived> 309 const typename eigen2_product_return_type<OtherDerived>::type 310 operator*(const EigenBase<OtherDerived>& rhs) const 311 { 312 typename OtherDerived::PlainObject::DenseType rhsPlainObject; 313 rhs.evalTo(rhsPlainObject); 314 return this->toDenseMatrix() * rhsPlainObject; 315 } 316 template<typename OtherMatrixType> 317 bool isApprox(const TriangularView<OtherMatrixType, Mode>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const 318 { 319 return this->toDenseMatrix().isApprox(other.toDenseMatrix(), precision); 320 } 321 template<typename OtherDerived> 322 bool isApprox(const MatrixBase<OtherDerived>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const 323 { 324 return this->toDenseMatrix().isApprox(other, precision); 325 } 326 #endif // EIGEN2_SUPPORT 327 328 template<int Side, typename Other> 329 inline const internal::triangular_solve_retval<Side,TriangularView, Other> 330 solve(const MatrixBase<Other>& other) const; 331 332 template<int Side, typename OtherDerived> 333 void solveInPlace(const MatrixBase<OtherDerived>& other) const; 334 335 template<typename Other> 336 inline const internal::triangular_solve_retval<OnTheLeft,TriangularView, Other> 337 solve(const MatrixBase<Other>& other) const 338 { return solve<OnTheLeft>(other); } 339 340 template<typename OtherDerived> 341 void solveInPlace(const MatrixBase<OtherDerived>& other) const 342 { return solveInPlace<OnTheLeft>(other); } 343 344 const SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() const 345 { 346 EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR); 347 return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix); 348 } 349 SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() 350 { 351 EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR); 352 return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix); 353 } 354 355 template<typename OtherDerived> 356 void swap(TriangularBase<OtherDerived> const & other) 357 { 358 TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived()); 359 } 360 361 template<typename OtherDerived> 362 void swap(MatrixBase<OtherDerived> const & other) 363 { 364 SwapWrapper<MatrixType> swaper(const_cast<MatrixType&>(m_matrix)); 365 TriangularView<SwapWrapper<MatrixType>,Mode>(swaper).lazyAssign(other.derived()); 366 } 367 368 Scalar determinant() const 369 { 370 if (Mode & UnitDiag) 371 return 1; 372 else if (Mode & ZeroDiag) 373 return 0; 374 else 375 return m_matrix.diagonal().prod(); 376 } 377 378 // TODO simplify the following: 379 template<typename ProductDerived, typename Lhs, typename Rhs> 380 EIGEN_STRONG_INLINE TriangularView& operator=(const ProductBase<ProductDerived, Lhs,Rhs>& other) 381 { 382 setZero(); 383 return assignProduct(other.derived(),1); 384 } 385 386 template<typename ProductDerived, typename Lhs, typename Rhs> 387 EIGEN_STRONG_INLINE TriangularView& operator+=(const ProductBase<ProductDerived, Lhs,Rhs>& other) 388 { 389 return assignProduct(other.derived(),1); 390 } 391 392 template<typename ProductDerived, typename Lhs, typename Rhs> 393 EIGEN_STRONG_INLINE TriangularView& operator-=(const ProductBase<ProductDerived, Lhs,Rhs>& other) 394 { 395 return assignProduct(other.derived(),-1); 396 } 397 398 399 template<typename ProductDerived> 400 EIGEN_STRONG_INLINE TriangularView& operator=(const ScaledProduct<ProductDerived>& other) 401 { 402 setZero(); 403 return assignProduct(other.derived(),other.alpha()); 404 } 405 406 template<typename ProductDerived> 407 EIGEN_STRONG_INLINE TriangularView& operator+=(const ScaledProduct<ProductDerived>& other) 408 { 409 return assignProduct(other.derived(),other.alpha()); 410 } 411 412 template<typename ProductDerived> 413 EIGEN_STRONG_INLINE TriangularView& operator-=(const ScaledProduct<ProductDerived>& other) 414 { 415 return assignProduct(other.derived(),-other.alpha()); 416 } 417 418 protected: 419 420 template<typename ProductDerived, typename Lhs, typename Rhs> 421 EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase<ProductDerived, Lhs,Rhs>& prod, const Scalar& alpha); 422 423 template<int Mode, bool LhsIsTriangular, 424 typename Lhs, bool LhsIsVector, 425 typename Rhs, bool RhsIsVector> 426 EIGEN_STRONG_INLINE TriangularView& assignProduct(const TriangularProduct<Mode, LhsIsTriangular, Lhs, LhsIsVector, Rhs, RhsIsVector>& prod, const Scalar& alpha) 427 { 428 lazyAssign(alpha*prod.eval()); 429 return *this; 430 } 431 432 MatrixTypeNested m_matrix; 433 }; 434 435 /*************************************************************************** 436 * Implementation of triangular evaluation/assignment 437 ***************************************************************************/ 438 439 namespace internal { 440 441 template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount, bool ClearOpposite> 442 struct triangular_assignment_selector 443 { 444 enum { 445 col = (UnrollCount-1) / Derived1::RowsAtCompileTime, 446 row = (UnrollCount-1) % Derived1::RowsAtCompileTime 447 }; 448 449 typedef typename Derived1::Scalar Scalar; 450 451 static inline void run(Derived1 &dst, const Derived2 &src) 452 { 453 triangular_assignment_selector<Derived1, Derived2, Mode, UnrollCount-1, ClearOpposite>::run(dst, src); 454 455 eigen_assert( Mode == Upper || Mode == Lower 456 || Mode == StrictlyUpper || Mode == StrictlyLower 457 || Mode == UnitUpper || Mode == UnitLower); 458 if((Mode == Upper && row <= col) 459 || (Mode == Lower && row >= col) 460 || (Mode == StrictlyUpper && row < col) 461 || (Mode == StrictlyLower && row > col) 462 || (Mode == UnitUpper && row < col) 463 || (Mode == UnitLower && row > col)) 464 dst.copyCoeff(row, col, src); 465 else if(ClearOpposite) 466 { 467 if (Mode&UnitDiag && row==col) 468 dst.coeffRef(row, col) = Scalar(1); 469 else 470 dst.coeffRef(row, col) = Scalar(0); 471 } 472 } 473 }; 474 475 // prevent buggy user code from causing an infinite recursion 476 template<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite> 477 struct triangular_assignment_selector<Derived1, Derived2, Mode, 0, ClearOpposite> 478 { 479 static inline void run(Derived1 &, const Derived2 &) {} 480 }; 481 482 template<typename Derived1, typename Derived2, bool ClearOpposite> 483 struct triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearOpposite> 484 { 485 typedef typename Derived1::Index Index; 486 typedef typename Derived1::Scalar Scalar; 487 static inline void run(Derived1 &dst, const Derived2 &src) 488 { 489 for(Index j = 0; j < dst.cols(); ++j) 490 { 491 Index maxi = (std::min)(j, dst.rows()-1); 492 for(Index i = 0; i <= maxi; ++i) 493 dst.copyCoeff(i, j, src); 494 if (ClearOpposite) 495 for(Index i = maxi+1; i < dst.rows(); ++i) 496 dst.coeffRef(i, j) = Scalar(0); 497 } 498 } 499 }; 500 501 template<typename Derived1, typename Derived2, bool ClearOpposite> 502 struct triangular_assignment_selector<Derived1, Derived2, Lower, Dynamic, ClearOpposite> 503 { 504 typedef typename Derived1::Index Index; 505 static inline void run(Derived1 &dst, const Derived2 &src) 506 { 507 for(Index j = 0; j < dst.cols(); ++j) 508 { 509 for(Index i = j; i < dst.rows(); ++i) 510 dst.copyCoeff(i, j, src); 511 Index maxi = (std::min)(j, dst.rows()); 512 if (ClearOpposite) 513 for(Index i = 0; i < maxi; ++i) 514 dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0); 515 } 516 } 517 }; 518 519 template<typename Derived1, typename Derived2, bool ClearOpposite> 520 struct triangular_assignment_selector<Derived1, Derived2, StrictlyUpper, Dynamic, ClearOpposite> 521 { 522 typedef typename Derived1::Index Index; 523 typedef typename Derived1::Scalar Scalar; 524 static inline void run(Derived1 &dst, const Derived2 &src) 525 { 526 for(Index j = 0; j < dst.cols(); ++j) 527 { 528 Index maxi = (std::min)(j, dst.rows()); 529 for(Index i = 0; i < maxi; ++i) 530 dst.copyCoeff(i, j, src); 531 if (ClearOpposite) 532 for(Index i = maxi; i < dst.rows(); ++i) 533 dst.coeffRef(i, j) = Scalar(0); 534 } 535 } 536 }; 537 538 template<typename Derived1, typename Derived2, bool ClearOpposite> 539 struct triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic, ClearOpposite> 540 { 541 typedef typename Derived1::Index Index; 542 static inline void run(Derived1 &dst, const Derived2 &src) 543 { 544 for(Index j = 0; j < dst.cols(); ++j) 545 { 546 for(Index i = j+1; i < dst.rows(); ++i) 547 dst.copyCoeff(i, j, src); 548 Index maxi = (std::min)(j, dst.rows()-1); 549 if (ClearOpposite) 550 for(Index i = 0; i <= maxi; ++i) 551 dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0); 552 } 553 } 554 }; 555 556 template<typename Derived1, typename Derived2, bool ClearOpposite> 557 struct triangular_assignment_selector<Derived1, Derived2, UnitUpper, Dynamic, ClearOpposite> 558 { 559 typedef typename Derived1::Index Index; 560 static inline void run(Derived1 &dst, const Derived2 &src) 561 { 562 for(Index j = 0; j < dst.cols(); ++j) 563 { 564 Index maxi = (std::min)(j, dst.rows()); 565 for(Index i = 0; i < maxi; ++i) 566 dst.copyCoeff(i, j, src); 567 if (ClearOpposite) 568 { 569 for(Index i = maxi+1; i < dst.rows(); ++i) 570 dst.coeffRef(i, j) = 0; 571 } 572 } 573 dst.diagonal().setOnes(); 574 } 575 }; 576 template<typename Derived1, typename Derived2, bool ClearOpposite> 577 struct triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, ClearOpposite> 578 { 579 typedef typename Derived1::Index Index; 580 static inline void run(Derived1 &dst, const Derived2 &src) 581 { 582 for(Index j = 0; j < dst.cols(); ++j) 583 { 584 Index maxi = (std::min)(j, dst.rows()); 585 for(Index i = maxi+1; i < dst.rows(); ++i) 586 dst.copyCoeff(i, j, src); 587 if (ClearOpposite) 588 { 589 for(Index i = 0; i < maxi; ++i) 590 dst.coeffRef(i, j) = 0; 591 } 592 } 593 dst.diagonal().setOnes(); 594 } 595 }; 596 597 } // end namespace internal 598 599 // FIXME should we keep that possibility 600 template<typename MatrixType, unsigned int Mode> 601 template<typename OtherDerived> 602 inline TriangularView<MatrixType, Mode>& 603 TriangularView<MatrixType, Mode>::operator=(const MatrixBase<OtherDerived>& other) 604 { 605 if(OtherDerived::Flags & EvalBeforeAssigningBit) 606 { 607 typename internal::plain_matrix_type<OtherDerived>::type other_evaluated(other.rows(), other.cols()); 608 other_evaluated.template triangularView<Mode>().lazyAssign(other.derived()); 609 lazyAssign(other_evaluated); 610 } 611 else 612 lazyAssign(other.derived()); 613 return *this; 614 } 615 616 // FIXME should we keep that possibility 617 template<typename MatrixType, unsigned int Mode> 618 template<typename OtherDerived> 619 void TriangularView<MatrixType, Mode>::lazyAssign(const MatrixBase<OtherDerived>& other) 620 { 621 enum { 622 unroll = MatrixType::SizeAtCompileTime != Dynamic 623 && internal::traits<OtherDerived>::CoeffReadCost != Dynamic 624 && MatrixType::SizeAtCompileTime*internal::traits<OtherDerived>::CoeffReadCost/2 <= EIGEN_UNROLLING_LIMIT 625 }; 626 eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols()); 627 628 internal::triangular_assignment_selector 629 <MatrixType, OtherDerived, int(Mode), 630 unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic, 631 false // do not change the opposite triangular part 632 >::run(m_matrix.const_cast_derived(), other.derived()); 633 } 634 635 636 637 template<typename MatrixType, unsigned int Mode> 638 template<typename OtherDerived> 639 inline TriangularView<MatrixType, Mode>& 640 TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>& other) 641 { 642 eigen_assert(Mode == int(OtherDerived::Mode)); 643 if(internal::traits<OtherDerived>::Flags & EvalBeforeAssigningBit) 644 { 645 typename OtherDerived::DenseMatrixType other_evaluated(other.rows(), other.cols()); 646 other_evaluated.template triangularView<Mode>().lazyAssign(other.derived().nestedExpression()); 647 lazyAssign(other_evaluated); 648 } 649 else 650 lazyAssign(other.derived().nestedExpression()); 651 return *this; 652 } 653 654 template<typename MatrixType, unsigned int Mode> 655 template<typename OtherDerived> 656 void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDerived>& other) 657 { 658 enum { 659 unroll = MatrixType::SizeAtCompileTime != Dynamic 660 && internal::traits<OtherDerived>::CoeffReadCost != Dynamic 661 && MatrixType::SizeAtCompileTime * internal::traits<OtherDerived>::CoeffReadCost / 2 662 <= EIGEN_UNROLLING_LIMIT 663 }; 664 eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols()); 665 666 internal::triangular_assignment_selector 667 <MatrixType, OtherDerived, int(Mode), 668 unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic, 669 false // preserve the opposite triangular part 670 >::run(m_matrix.const_cast_derived(), other.derived().nestedExpression()); 671 } 672 673 /*************************************************************************** 674 * Implementation of TriangularBase methods 675 ***************************************************************************/ 676 677 /** Assigns a triangular or selfadjoint matrix to a dense matrix. 678 * If the matrix is triangular, the opposite part is set to zero. */ 679 template<typename Derived> 680 template<typename DenseDerived> 681 void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const 682 { 683 if(internal::traits<Derived>::Flags & EvalBeforeAssigningBit) 684 { 685 typename internal::plain_matrix_type<Derived>::type other_evaluated(rows(), cols()); 686 evalToLazy(other_evaluated); 687 other.derived().swap(other_evaluated); 688 } 689 else 690 evalToLazy(other.derived()); 691 } 692 693 /** Assigns a triangular or selfadjoint matrix to a dense matrix. 694 * If the matrix is triangular, the opposite part is set to zero. */ 695 template<typename Derived> 696 template<typename DenseDerived> 697 void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const 698 { 699 enum { 700 unroll = DenseDerived::SizeAtCompileTime != Dynamic 701 && internal::traits<Derived>::CoeffReadCost != Dynamic 702 && DenseDerived::SizeAtCompileTime * internal::traits<Derived>::CoeffReadCost / 2 703 <= EIGEN_UNROLLING_LIMIT 704 }; 705 other.derived().resize(this->rows(), this->cols()); 706 707 internal::triangular_assignment_selector 708 <DenseDerived, typename internal::traits<Derived>::MatrixTypeNestedCleaned, Derived::Mode, 709 unroll ? int(DenseDerived::SizeAtCompileTime) : Dynamic, 710 true // clear the opposite triangular part 711 >::run(other.derived(), derived().nestedExpression()); 712 } 713 714 /*************************************************************************** 715 * Implementation of TriangularView methods 716 ***************************************************************************/ 717 718 /*************************************************************************** 719 * Implementation of MatrixBase methods 720 ***************************************************************************/ 721 722 #ifdef EIGEN2_SUPPORT 723 724 // implementation of part<>(), including the SelfAdjoint case. 725 726 namespace internal { 727 template<typename MatrixType, unsigned int Mode> 728 struct eigen2_part_return_type 729 { 730 typedef TriangularView<MatrixType, Mode> type; 731 }; 732 733 template<typename MatrixType> 734 struct eigen2_part_return_type<MatrixType, SelfAdjoint> 735 { 736 typedef SelfAdjointView<MatrixType, Upper> type; 737 }; 738 } 739 740 /** \deprecated use MatrixBase::triangularView() */ 741 template<typename Derived> 742 template<unsigned int Mode> 743 const typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part() const 744 { 745 return derived(); 746 } 747 748 /** \deprecated use MatrixBase::triangularView() */ 749 template<typename Derived> 750 template<unsigned int Mode> 751 typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part() 752 { 753 return derived(); 754 } 755 #endif 756 757 /** 758 * \returns an expression of a triangular view extracted from the current matrix 759 * 760 * The parameter \a Mode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper, 761 * \c #Lower, \c #StrictlyLower, \c #UnitLower. 762 * 763 * Example: \include MatrixBase_extract.cpp 764 * Output: \verbinclude MatrixBase_extract.out 765 * 766 * \sa class TriangularView 767 */ 768 template<typename Derived> 769 template<unsigned int Mode> 770 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type 771 MatrixBase<Derived>::triangularView() 772 { 773 return derived(); 774 } 775 776 /** This is the const version of MatrixBase::triangularView() */ 777 template<typename Derived> 778 template<unsigned int Mode> 779 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type 780 MatrixBase<Derived>::triangularView() const 781 { 782 return derived(); 783 } 784 785 /** \returns true if *this is approximately equal to an upper triangular matrix, 786 * within the precision given by \a prec. 787 * 788 * \sa isLowerTriangular() 789 */ 790 template<typename Derived> 791 bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const 792 { 793 using std::abs; 794 RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1); 795 for(Index j = 0; j < cols(); ++j) 796 { 797 Index maxi = (std::min)(j, rows()-1); 798 for(Index i = 0; i <= maxi; ++i) 799 { 800 RealScalar absValue = abs(coeff(i,j)); 801 if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue; 802 } 803 } 804 RealScalar threshold = maxAbsOnUpperPart * prec; 805 for(Index j = 0; j < cols(); ++j) 806 for(Index i = j+1; i < rows(); ++i) 807 if(abs(coeff(i, j)) > threshold) return false; 808 return true; 809 } 810 811 /** \returns true if *this is approximately equal to a lower triangular matrix, 812 * within the precision given by \a prec. 813 * 814 * \sa isUpperTriangular() 815 */ 816 template<typename Derived> 817 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const 818 { 819 using std::abs; 820 RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1); 821 for(Index j = 0; j < cols(); ++j) 822 for(Index i = j; i < rows(); ++i) 823 { 824 RealScalar absValue = abs(coeff(i,j)); 825 if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue; 826 } 827 RealScalar threshold = maxAbsOnLowerPart * prec; 828 for(Index j = 1; j < cols(); ++j) 829 { 830 Index maxi = (std::min)(j, rows()-1); 831 for(Index i = 0; i < maxi; ++i) 832 if(abs(coeff(i, j)) > threshold) return false; 833 } 834 return true; 835 } 836 837 } // end namespace Eigen 838 839 #endif // EIGEN_TRIANGULARMATRIX_H 840
