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,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,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,-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,other.alpha()); 404 } 405 406 template<typename ProductDerived> 407 EIGEN_STRONG_INLINE TriangularView& operator+=(const ScaledProduct<ProductDerived>& other) 408 { 409 return assignProduct(other,other.alpha()); 410 } 411 412 template<typename ProductDerived> 413 EIGEN_STRONG_INLINE TriangularView& operator-=(const ScaledProduct<ProductDerived>& other) 414 { 415 return assignProduct(other,-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 MatrixTypeNested m_matrix; 424 }; 425 426 /*************************************************************************** 427 * Implementation of triangular evaluation/assignment 428 ***************************************************************************/ 429 430 namespace internal { 431 432 template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount, bool ClearOpposite> 433 struct triangular_assignment_selector 434 { 435 enum { 436 col = (UnrollCount-1) / Derived1::RowsAtCompileTime, 437 row = (UnrollCount-1) % Derived1::RowsAtCompileTime 438 }; 439 440 typedef typename Derived1::Scalar Scalar; 441 442 static inline void run(Derived1 &dst, const Derived2 &src) 443 { 444 triangular_assignment_selector<Derived1, Derived2, Mode, UnrollCount-1, ClearOpposite>::run(dst, src); 445 446 eigen_assert( Mode == Upper || Mode == Lower 447 || Mode == StrictlyUpper || Mode == StrictlyLower 448 || Mode == UnitUpper || Mode == UnitLower); 449 if((Mode == Upper && row <= col) 450 || (Mode == Lower && row >= col) 451 || (Mode == StrictlyUpper && row < col) 452 || (Mode == StrictlyLower && row > col) 453 || (Mode == UnitUpper && row < col) 454 || (Mode == UnitLower && row > col)) 455 dst.copyCoeff(row, col, src); 456 else if(ClearOpposite) 457 { 458 if (Mode&UnitDiag && row==col) 459 dst.coeffRef(row, col) = Scalar(1); 460 else 461 dst.coeffRef(row, col) = Scalar(0); 462 } 463 } 464 }; 465 466 // prevent buggy user code from causing an infinite recursion 467 template<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite> 468 struct triangular_assignment_selector<Derived1, Derived2, Mode, 0, ClearOpposite> 469 { 470 static inline void run(Derived1 &, const Derived2 &) {} 471 }; 472 473 template<typename Derived1, typename Derived2, bool ClearOpposite> 474 struct triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearOpposite> 475 { 476 typedef typename Derived1::Index Index; 477 typedef typename Derived1::Scalar Scalar; 478 static inline void run(Derived1 &dst, const Derived2 &src) 479 { 480 for(Index j = 0; j < dst.cols(); ++j) 481 { 482 Index maxi = (std::min)(j, dst.rows()-1); 483 for(Index i = 0; i <= maxi; ++i) 484 dst.copyCoeff(i, j, src); 485 if (ClearOpposite) 486 for(Index i = maxi+1; i < dst.rows(); ++i) 487 dst.coeffRef(i, j) = Scalar(0); 488 } 489 } 490 }; 491 492 template<typename Derived1, typename Derived2, bool ClearOpposite> 493 struct triangular_assignment_selector<Derived1, Derived2, Lower, Dynamic, ClearOpposite> 494 { 495 typedef typename Derived1::Index Index; 496 static inline void run(Derived1 &dst, const Derived2 &src) 497 { 498 for(Index j = 0; j < dst.cols(); ++j) 499 { 500 for(Index i = j; i < dst.rows(); ++i) 501 dst.copyCoeff(i, j, src); 502 Index maxi = (std::min)(j, dst.rows()); 503 if (ClearOpposite) 504 for(Index i = 0; i < maxi; ++i) 505 dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0); 506 } 507 } 508 }; 509 510 template<typename Derived1, typename Derived2, bool ClearOpposite> 511 struct triangular_assignment_selector<Derived1, Derived2, StrictlyUpper, Dynamic, ClearOpposite> 512 { 513 typedef typename Derived1::Index Index; 514 typedef typename Derived1::Scalar Scalar; 515 static inline void run(Derived1 &dst, const Derived2 &src) 516 { 517 for(Index j = 0; j < dst.cols(); ++j) 518 { 519 Index maxi = (std::min)(j, dst.rows()); 520 for(Index i = 0; i < maxi; ++i) 521 dst.copyCoeff(i, j, src); 522 if (ClearOpposite) 523 for(Index i = maxi; i < dst.rows(); ++i) 524 dst.coeffRef(i, j) = Scalar(0); 525 } 526 } 527 }; 528 529 template<typename Derived1, typename Derived2, bool ClearOpposite> 530 struct triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic, ClearOpposite> 531 { 532 typedef typename Derived1::Index Index; 533 static inline void run(Derived1 &dst, const Derived2 &src) 534 { 535 for(Index j = 0; j < dst.cols(); ++j) 536 { 537 for(Index i = j+1; i < dst.rows(); ++i) 538 dst.copyCoeff(i, j, src); 539 Index maxi = (std::min)(j, dst.rows()-1); 540 if (ClearOpposite) 541 for(Index i = 0; i <= maxi; ++i) 542 dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0); 543 } 544 } 545 }; 546 547 template<typename Derived1, typename Derived2, bool ClearOpposite> 548 struct triangular_assignment_selector<Derived1, Derived2, UnitUpper, Dynamic, ClearOpposite> 549 { 550 typedef typename Derived1::Index Index; 551 static inline void run(Derived1 &dst, const Derived2 &src) 552 { 553 for(Index j = 0; j < dst.cols(); ++j) 554 { 555 Index maxi = (std::min)(j, dst.rows()); 556 for(Index i = 0; i < maxi; ++i) 557 dst.copyCoeff(i, j, src); 558 if (ClearOpposite) 559 { 560 for(Index i = maxi+1; i < dst.rows(); ++i) 561 dst.coeffRef(i, j) = 0; 562 } 563 } 564 dst.diagonal().setOnes(); 565 } 566 }; 567 template<typename Derived1, typename Derived2, bool ClearOpposite> 568 struct triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, ClearOpposite> 569 { 570 typedef typename Derived1::Index Index; 571 static inline void run(Derived1 &dst, const Derived2 &src) 572 { 573 for(Index j = 0; j < dst.cols(); ++j) 574 { 575 Index maxi = (std::min)(j, dst.rows()); 576 for(Index i = maxi+1; i < dst.rows(); ++i) 577 dst.copyCoeff(i, j, src); 578 if (ClearOpposite) 579 { 580 for(Index i = 0; i < maxi; ++i) 581 dst.coeffRef(i, j) = 0; 582 } 583 } 584 dst.diagonal().setOnes(); 585 } 586 }; 587 588 } // end namespace internal 589 590 // FIXME should we keep that possibility 591 template<typename MatrixType, unsigned int Mode> 592 template<typename OtherDerived> 593 inline TriangularView<MatrixType, Mode>& 594 TriangularView<MatrixType, Mode>::operator=(const MatrixBase<OtherDerived>& other) 595 { 596 if(OtherDerived::Flags & EvalBeforeAssigningBit) 597 { 598 typename internal::plain_matrix_type<OtherDerived>::type other_evaluated(other.rows(), other.cols()); 599 other_evaluated.template triangularView<Mode>().lazyAssign(other.derived()); 600 lazyAssign(other_evaluated); 601 } 602 else 603 lazyAssign(other.derived()); 604 return *this; 605 } 606 607 // FIXME should we keep that possibility 608 template<typename MatrixType, unsigned int Mode> 609 template<typename OtherDerived> 610 void TriangularView<MatrixType, Mode>::lazyAssign(const MatrixBase<OtherDerived>& other) 611 { 612 enum { 613 unroll = MatrixType::SizeAtCompileTime != Dynamic 614 && internal::traits<OtherDerived>::CoeffReadCost != Dynamic 615 && MatrixType::SizeAtCompileTime*internal::traits<OtherDerived>::CoeffReadCost/2 <= EIGEN_UNROLLING_LIMIT 616 }; 617 eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols()); 618 619 internal::triangular_assignment_selector 620 <MatrixType, OtherDerived, int(Mode), 621 unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic, 622 false // do not change the opposite triangular part 623 >::run(m_matrix.const_cast_derived(), other.derived()); 624 } 625 626 627 628 template<typename MatrixType, unsigned int Mode> 629 template<typename OtherDerived> 630 inline TriangularView<MatrixType, Mode>& 631 TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>& other) 632 { 633 eigen_assert(Mode == int(OtherDerived::Mode)); 634 if(internal::traits<OtherDerived>::Flags & EvalBeforeAssigningBit) 635 { 636 typename OtherDerived::DenseMatrixType other_evaluated(other.rows(), other.cols()); 637 other_evaluated.template triangularView<Mode>().lazyAssign(other.derived().nestedExpression()); 638 lazyAssign(other_evaluated); 639 } 640 else 641 lazyAssign(other.derived().nestedExpression()); 642 return *this; 643 } 644 645 template<typename MatrixType, unsigned int Mode> 646 template<typename OtherDerived> 647 void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDerived>& other) 648 { 649 enum { 650 unroll = MatrixType::SizeAtCompileTime != Dynamic 651 && internal::traits<OtherDerived>::CoeffReadCost != Dynamic 652 && MatrixType::SizeAtCompileTime * internal::traits<OtherDerived>::CoeffReadCost / 2 653 <= EIGEN_UNROLLING_LIMIT 654 }; 655 eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols()); 656 657 internal::triangular_assignment_selector 658 <MatrixType, OtherDerived, int(Mode), 659 unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic, 660 false // preserve the opposite triangular part 661 >::run(m_matrix.const_cast_derived(), other.derived().nestedExpression()); 662 } 663 664 /*************************************************************************** 665 * Implementation of TriangularBase methods 666 ***************************************************************************/ 667 668 /** Assigns a triangular or selfadjoint matrix to a dense matrix. 669 * If the matrix is triangular, the opposite part is set to zero. */ 670 template<typename Derived> 671 template<typename DenseDerived> 672 void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const 673 { 674 if(internal::traits<Derived>::Flags & EvalBeforeAssigningBit) 675 { 676 typename internal::plain_matrix_type<Derived>::type other_evaluated(rows(), cols()); 677 evalToLazy(other_evaluated); 678 other.derived().swap(other_evaluated); 679 } 680 else 681 evalToLazy(other.derived()); 682 } 683 684 /** Assigns a triangular or selfadjoint matrix to a dense matrix. 685 * If the matrix is triangular, the opposite part is set to zero. */ 686 template<typename Derived> 687 template<typename DenseDerived> 688 void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const 689 { 690 enum { 691 unroll = DenseDerived::SizeAtCompileTime != Dynamic 692 && internal::traits<Derived>::CoeffReadCost != Dynamic 693 && DenseDerived::SizeAtCompileTime * internal::traits<Derived>::CoeffReadCost / 2 694 <= EIGEN_UNROLLING_LIMIT 695 }; 696 other.derived().resize(this->rows(), this->cols()); 697 698 internal::triangular_assignment_selector 699 <DenseDerived, typename internal::traits<Derived>::MatrixTypeNestedCleaned, Derived::Mode, 700 unroll ? int(DenseDerived::SizeAtCompileTime) : Dynamic, 701 true // clear the opposite triangular part 702 >::run(other.derived(), derived().nestedExpression()); 703 } 704 705 /*************************************************************************** 706 * Implementation of TriangularView methods 707 ***************************************************************************/ 708 709 /*************************************************************************** 710 * Implementation of MatrixBase methods 711 ***************************************************************************/ 712 713 #ifdef EIGEN2_SUPPORT 714 715 // implementation of part<>(), including the SelfAdjoint case. 716 717 namespace internal { 718 template<typename MatrixType, unsigned int Mode> 719 struct eigen2_part_return_type 720 { 721 typedef TriangularView<MatrixType, Mode> type; 722 }; 723 724 template<typename MatrixType> 725 struct eigen2_part_return_type<MatrixType, SelfAdjoint> 726 { 727 typedef SelfAdjointView<MatrixType, Upper> type; 728 }; 729 } 730 731 /** \deprecated use MatrixBase::triangularView() */ 732 template<typename Derived> 733 template<unsigned int Mode> 734 const typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part() const 735 { 736 return derived(); 737 } 738 739 /** \deprecated use MatrixBase::triangularView() */ 740 template<typename Derived> 741 template<unsigned int Mode> 742 typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part() 743 { 744 return derived(); 745 } 746 #endif 747 748 /** 749 * \returns an expression of a triangular view extracted from the current matrix 750 * 751 * The parameter \a Mode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper, 752 * \c #Lower, \c #StrictlyLower, \c #UnitLower. 753 * 754 * Example: \include MatrixBase_extract.cpp 755 * Output: \verbinclude MatrixBase_extract.out 756 * 757 * \sa class TriangularView 758 */ 759 template<typename Derived> 760 template<unsigned int Mode> 761 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type 762 MatrixBase<Derived>::triangularView() 763 { 764 return derived(); 765 } 766 767 /** This is the const version of MatrixBase::triangularView() */ 768 template<typename Derived> 769 template<unsigned int Mode> 770 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type 771 MatrixBase<Derived>::triangularView() const 772 { 773 return derived(); 774 } 775 776 /** \returns true if *this is approximately equal to an upper triangular matrix, 777 * within the precision given by \a prec. 778 * 779 * \sa isLowerTriangular() 780 */ 781 template<typename Derived> 782 bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const 783 { 784 using std::abs; 785 RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1); 786 for(Index j = 0; j < cols(); ++j) 787 { 788 Index maxi = (std::min)(j, rows()-1); 789 for(Index i = 0; i <= maxi; ++i) 790 { 791 RealScalar absValue = abs(coeff(i,j)); 792 if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue; 793 } 794 } 795 RealScalar threshold = maxAbsOnUpperPart * prec; 796 for(Index j = 0; j < cols(); ++j) 797 for(Index i = j+1; i < rows(); ++i) 798 if(abs(coeff(i, j)) > threshold) return false; 799 return true; 800 } 801 802 /** \returns true if *this is approximately equal to a lower triangular matrix, 803 * within the precision given by \a prec. 804 * 805 * \sa isUpperTriangular() 806 */ 807 template<typename Derived> 808 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const 809 { 810 using std::abs; 811 RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1); 812 for(Index j = 0; j < cols(); ++j) 813 for(Index i = j; i < rows(); ++i) 814 { 815 RealScalar absValue = abs(coeff(i,j)); 816 if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue; 817 } 818 RealScalar threshold = maxAbsOnLowerPart * prec; 819 for(Index j = 1; j < cols(); ++j) 820 { 821 Index maxi = (std::min)(j, rows()-1); 822 for(Index i = 0; i < maxi; ++i) 823 if(abs(coeff(i, j)) > threshold) return false; 824 } 825 return true; 826 } 827 828 } // end namespace Eigen 829 830 #endif // EIGEN_TRIANGULARMATRIX_H 831