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@gmail.com> 5 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@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 TriangularBase()47 inline TriangularBase() { eigen_assert(!((Mode&UnitDiag) && (Mode&ZeroDiag))); } 48 rows()49 inline Index rows() const { return derived().rows(); } cols()50 inline Index cols() const { return derived().cols(); } outerStride()51 inline Index outerStride() const { return derived().outerStride(); } innerStride()52 inline Index innerStride() const { return derived().innerStride(); } 53 coeff(Index row,Index col)54 inline Scalar coeff(Index row, Index col) const { return derived().coeff(row,col); } coeffRef(Index row,Index 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> copyCoeff(Index row,Index col,Other & 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 operator()65 inline Scalar operator()(Index row, Index col) const 66 { 67 check_coordinates(row, col); 68 return coeff(row,col); 69 } operator()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 derived()77 inline const Derived& derived() const { return *static_cast<const Derived*>(this); } derived()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 toDenseMatrix()86 DenseMatrixType toDenseMatrix() const 87 { 88 DenseMatrixType res(rows(), cols()); 89 evalToLazy(res); 90 return res; 91 } 92 93 protected: 94 check_coordinates(Index row,Index col)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 check_coordinates_internal(Index row,Index col)109 void check_coordinates_internal(Index row, Index col) const 110 { 111 check_coordinates(row, col); 112 } 113 #else check_coordinates_internal(Index,Index)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::IsVectorAtCompileTime> 282 operator*(const MatrixBase<OtherDerived>& rhs) const 283 { 284 return TriangularProduct 285 <Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime> 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::IsVectorAtCompileTime,MatrixType,false> 292 operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs) 293 { 294 return TriangularProduct 295 <Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,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 static inline void run(Derived1 &dst, const Derived2 &src) 515 { 516 for(Index j = 0; j < dst.cols(); ++j) 517 { 518 Index maxi = (std::min)(j, dst.rows()); 519 for(Index i = 0; i < maxi; ++i) 520 dst.copyCoeff(i, j, src); 521 if (ClearOpposite) 522 for(Index i = maxi; i < dst.rows(); ++i) 523 dst.coeffRef(i, j) = 0; 524 } 525 } 526 }; 527 528 template<typename Derived1, typename Derived2, bool ClearOpposite> 529 struct triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic, ClearOpposite> 530 { 531 typedef typename Derived1::Index Index; 532 static inline void run(Derived1 &dst, const Derived2 &src) 533 { 534 for(Index j = 0; j < dst.cols(); ++j) 535 { 536 for(Index i = j+1; i < dst.rows(); ++i) 537 dst.copyCoeff(i, j, src); 538 Index maxi = (std::min)(j, dst.rows()-1); 539 if (ClearOpposite) 540 for(Index i = 0; i <= maxi; ++i) 541 dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0); 542 } 543 } 544 }; 545 546 template<typename Derived1, typename Derived2, bool ClearOpposite> 547 struct triangular_assignment_selector<Derived1, Derived2, UnitUpper, Dynamic, ClearOpposite> 548 { 549 typedef typename Derived1::Index Index; 550 static inline void run(Derived1 &dst, const Derived2 &src) 551 { 552 for(Index j = 0; j < dst.cols(); ++j) 553 { 554 Index maxi = (std::min)(j, dst.rows()); 555 for(Index i = 0; i < maxi; ++i) 556 dst.copyCoeff(i, j, src); 557 if (ClearOpposite) 558 { 559 for(Index i = maxi+1; i < dst.rows(); ++i) 560 dst.coeffRef(i, j) = 0; 561 } 562 } 563 dst.diagonal().setOnes(); 564 } 565 }; 566 template<typename Derived1, typename Derived2, bool ClearOpposite> 567 struct triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, ClearOpposite> 568 { 569 typedef typename Derived1::Index Index; 570 static inline void run(Derived1 &dst, const Derived2 &src) 571 { 572 for(Index j = 0; j < dst.cols(); ++j) 573 { 574 Index maxi = (std::min)(j, dst.rows()); 575 for(Index i = maxi+1; i < dst.rows(); ++i) 576 dst.copyCoeff(i, j, src); 577 if (ClearOpposite) 578 { 579 for(Index i = 0; i < maxi; ++i) 580 dst.coeffRef(i, j) = 0; 581 } 582 } 583 dst.diagonal().setOnes(); 584 } 585 }; 586 587 } // end namespace internal 588 589 // FIXME should we keep that possibility 590 template<typename MatrixType, unsigned int Mode> 591 template<typename OtherDerived> 592 inline TriangularView<MatrixType, Mode>& 593 TriangularView<MatrixType, Mode>::operator=(const MatrixBase<OtherDerived>& other) 594 { 595 if(OtherDerived::Flags & EvalBeforeAssigningBit) 596 { 597 typename internal::plain_matrix_type<OtherDerived>::type other_evaluated(other.rows(), other.cols()); 598 other_evaluated.template triangularView<Mode>().lazyAssign(other.derived()); 599 lazyAssign(other_evaluated); 600 } 601 else 602 lazyAssign(other.derived()); 603 return *this; 604 } 605 606 // FIXME should we keep that possibility 607 template<typename MatrixType, unsigned int Mode> 608 template<typename OtherDerived> 609 void TriangularView<MatrixType, Mode>::lazyAssign(const MatrixBase<OtherDerived>& other) 610 { 611 enum { 612 unroll = MatrixType::SizeAtCompileTime != Dynamic 613 && internal::traits<OtherDerived>::CoeffReadCost != Dynamic 614 && MatrixType::SizeAtCompileTime*internal::traits<OtherDerived>::CoeffReadCost/2 <= EIGEN_UNROLLING_LIMIT 615 }; 616 eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols()); 617 618 internal::triangular_assignment_selector 619 <MatrixType, OtherDerived, int(Mode), 620 unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic, 621 false // do not change the opposite triangular part 622 >::run(m_matrix.const_cast_derived(), other.derived()); 623 } 624 625 626 627 template<typename MatrixType, unsigned int Mode> 628 template<typename OtherDerived> 629 inline TriangularView<MatrixType, Mode>& 630 TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>& other) 631 { 632 eigen_assert(Mode == int(OtherDerived::Mode)); 633 if(internal::traits<OtherDerived>::Flags & EvalBeforeAssigningBit) 634 { 635 typename OtherDerived::DenseMatrixType other_evaluated(other.rows(), other.cols()); 636 other_evaluated.template triangularView<Mode>().lazyAssign(other.derived().nestedExpression()); 637 lazyAssign(other_evaluated); 638 } 639 else 640 lazyAssign(other.derived().nestedExpression()); 641 return *this; 642 } 643 644 template<typename MatrixType, unsigned int Mode> 645 template<typename OtherDerived> 646 void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDerived>& other) 647 { 648 enum { 649 unroll = MatrixType::SizeAtCompileTime != Dynamic 650 && internal::traits<OtherDerived>::CoeffReadCost != Dynamic 651 && MatrixType::SizeAtCompileTime * internal::traits<OtherDerived>::CoeffReadCost / 2 652 <= EIGEN_UNROLLING_LIMIT 653 }; 654 eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols()); 655 656 internal::triangular_assignment_selector 657 <MatrixType, OtherDerived, int(Mode), 658 unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic, 659 false // preserve the opposite triangular part 660 >::run(m_matrix.const_cast_derived(), other.derived().nestedExpression()); 661 } 662 663 /*************************************************************************** 664 * Implementation of TriangularBase methods 665 ***************************************************************************/ 666 667 /** Assigns a triangular or selfadjoint matrix to a dense matrix. 668 * If the matrix is triangular, the opposite part is set to zero. */ 669 template<typename Derived> 670 template<typename DenseDerived> 671 void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const 672 { 673 if(internal::traits<Derived>::Flags & EvalBeforeAssigningBit) 674 { 675 typename internal::plain_matrix_type<Derived>::type other_evaluated(rows(), cols()); 676 evalToLazy(other_evaluated); 677 other.derived().swap(other_evaluated); 678 } 679 else 680 evalToLazy(other.derived()); 681 } 682 683 /** Assigns a triangular or selfadjoint matrix to a dense matrix. 684 * If the matrix is triangular, the opposite part is set to zero. */ 685 template<typename Derived> 686 template<typename DenseDerived> 687 void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const 688 { 689 enum { 690 unroll = DenseDerived::SizeAtCompileTime != Dynamic 691 && internal::traits<Derived>::CoeffReadCost != Dynamic 692 && DenseDerived::SizeAtCompileTime * internal::traits<Derived>::CoeffReadCost / 2 693 <= EIGEN_UNROLLING_LIMIT 694 }; 695 other.derived().resize(this->rows(), this->cols()); 696 697 internal::triangular_assignment_selector 698 <DenseDerived, typename internal::traits<Derived>::MatrixTypeNestedCleaned, Derived::Mode, 699 unroll ? int(DenseDerived::SizeAtCompileTime) : Dynamic, 700 true // clear the opposite triangular part 701 >::run(other.derived(), derived().nestedExpression()); 702 } 703 704 /*************************************************************************** 705 * Implementation of TriangularView methods 706 ***************************************************************************/ 707 708 /*************************************************************************** 709 * Implementation of MatrixBase methods 710 ***************************************************************************/ 711 712 #ifdef EIGEN2_SUPPORT 713 714 // implementation of part<>(), including the SelfAdjoint case. 715 716 namespace internal { 717 template<typename MatrixType, unsigned int Mode> 718 struct eigen2_part_return_type 719 { 720 typedef TriangularView<MatrixType, Mode> type; 721 }; 722 723 template<typename MatrixType> 724 struct eigen2_part_return_type<MatrixType, SelfAdjoint> 725 { 726 typedef SelfAdjointView<MatrixType, Upper> type; 727 }; 728 } 729 730 /** \deprecated use MatrixBase::triangularView() */ 731 template<typename Derived> 732 template<unsigned int Mode> 733 const typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part() const 734 { 735 return derived(); 736 } 737 738 /** \deprecated use MatrixBase::triangularView() */ 739 template<typename Derived> 740 template<unsigned int Mode> 741 typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part() 742 { 743 return derived(); 744 } 745 #endif 746 747 /** 748 * \returns an expression of a triangular view extracted from the current matrix 749 * 750 * The parameter \a Mode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper, 751 * \c #Lower, \c #StrictlyLower, \c #UnitLower. 752 * 753 * Example: \include MatrixBase_extract.cpp 754 * Output: \verbinclude MatrixBase_extract.out 755 * 756 * \sa class TriangularView 757 */ 758 template<typename Derived> 759 template<unsigned int Mode> 760 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type 761 MatrixBase<Derived>::triangularView() 762 { 763 return derived(); 764 } 765 766 /** This is the const version of MatrixBase::triangularView() */ 767 template<typename Derived> 768 template<unsigned int Mode> 769 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type 770 MatrixBase<Derived>::triangularView() const 771 { 772 return derived(); 773 } 774 775 /** \returns true if *this is approximately equal to an upper triangular matrix, 776 * within the precision given by \a prec. 777 * 778 * \sa isLowerTriangular() 779 */ 780 template<typename Derived> 781 bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const 782 { 783 RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1); 784 for(Index j = 0; j < cols(); ++j) 785 { 786 Index maxi = (std::min)(j, rows()-1); 787 for(Index i = 0; i <= maxi; ++i) 788 { 789 RealScalar absValue = internal::abs(coeff(i,j)); 790 if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue; 791 } 792 } 793 RealScalar threshold = maxAbsOnUpperPart * prec; 794 for(Index j = 0; j < cols(); ++j) 795 for(Index i = j+1; i < rows(); ++i) 796 if(internal::abs(coeff(i, j)) > threshold) return false; 797 return true; 798 } 799 800 /** \returns true if *this is approximately equal to a lower triangular matrix, 801 * within the precision given by \a prec. 802 * 803 * \sa isUpperTriangular() 804 */ 805 template<typename Derived> 806 bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const 807 { 808 RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1); 809 for(Index j = 0; j < cols(); ++j) 810 for(Index i = j; i < rows(); ++i) 811 { 812 RealScalar absValue = internal::abs(coeff(i,j)); 813 if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue; 814 } 815 RealScalar threshold = maxAbsOnLowerPart * prec; 816 for(Index j = 1; j < cols(); ++j) 817 { 818 Index maxi = (std::min)(j, rows()-1); 819 for(Index i = 0; i < maxi; ++i) 820 if(internal::abs(coeff(i, j)) > threshold) return false; 821 } 822 return true; 823 } 824 825 } // end namespace Eigen 826 827 #endif // EIGEN_TRIANGULARMATRIX_H 828