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 /** \class TriangularBase 23 * \ingroup Core_Module 24 * 25 * \brief Base class for triangular part in a matrix 26 */ 27 template<typename Derived> class TriangularBase : public EigenBase<Derived> 28 { 29 public: 30 31 enum { 32 Mode = internal::traits<Derived>::Mode, 33 RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime, 34 ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime, 35 MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime, 36 MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime, 37 38 SizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::RowsAtCompileTime, 39 internal::traits<Derived>::ColsAtCompileTime>::ret), 40 /**< This is equal to the number of coefficients, i.e. the number of 41 * rows times the number of columns, or to \a Dynamic if this is not 42 * known at compile-time. \sa RowsAtCompileTime, ColsAtCompileTime */ 43 44 MaxSizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::MaxRowsAtCompileTime, 45 internal::traits<Derived>::MaxColsAtCompileTime>::ret) 46 47 }; 48 typedef typename internal::traits<Derived>::Scalar Scalar; 49 typedef typename internal::traits<Derived>::StorageKind StorageKind; 50 typedef typename internal::traits<Derived>::StorageIndex StorageIndex; 51 typedef typename internal::traits<Derived>::FullMatrixType DenseMatrixType; 52 typedef DenseMatrixType DenseType; 53 typedef Derived const& Nested; 54 55 EIGEN_DEVICE_FUNC 56 inline TriangularBase() { eigen_assert(!((Mode&UnitDiag) && (Mode&ZeroDiag))); } 57 58 EIGEN_DEVICE_FUNC 59 inline Index rows() const { return derived().rows(); } 60 EIGEN_DEVICE_FUNC 61 inline Index cols() const { return derived().cols(); } 62 EIGEN_DEVICE_FUNC 63 inline Index outerStride() const { return derived().outerStride(); } 64 EIGEN_DEVICE_FUNC 65 inline Index innerStride() const { return derived().innerStride(); } 66 67 // dummy resize function 68 void resize(Index rows, Index cols) 69 { 70 EIGEN_UNUSED_VARIABLE(rows); 71 EIGEN_UNUSED_VARIABLE(cols); 72 eigen_assert(rows==this->rows() && cols==this->cols()); 73 } 74 75 EIGEN_DEVICE_FUNC 76 inline Scalar coeff(Index row, Index col) const { return derived().coeff(row,col); } 77 EIGEN_DEVICE_FUNC 78 inline Scalar& coeffRef(Index row, Index col) { return derived().coeffRef(row,col); } 79 80 /** \see MatrixBase::copyCoeff(row,col) 81 */ 82 template<typename Other> 83 EIGEN_DEVICE_FUNC 84 EIGEN_STRONG_INLINE void copyCoeff(Index row, Index col, Other& other) 85 { 86 derived().coeffRef(row, col) = other.coeff(row, col); 87 } 88 89 EIGEN_DEVICE_FUNC 90 inline Scalar operator()(Index row, Index col) const 91 { 92 check_coordinates(row, col); 93 return coeff(row,col); 94 } 95 EIGEN_DEVICE_FUNC 96 inline Scalar& operator()(Index row, Index col) 97 { 98 check_coordinates(row, col); 99 return coeffRef(row,col); 100 } 101 102 #ifndef EIGEN_PARSED_BY_DOXYGEN 103 EIGEN_DEVICE_FUNC 104 inline const Derived& derived() const { return *static_cast<const Derived*>(this); } 105 EIGEN_DEVICE_FUNC 106 inline Derived& derived() { return *static_cast<Derived*>(this); } 107 #endif // not EIGEN_PARSED_BY_DOXYGEN 108 109 template<typename DenseDerived> 110 EIGEN_DEVICE_FUNC 111 void evalTo(MatrixBase<DenseDerived> &other) const; 112 template<typename DenseDerived> 113 EIGEN_DEVICE_FUNC 114 void evalToLazy(MatrixBase<DenseDerived> &other) const; 115 116 EIGEN_DEVICE_FUNC 117 DenseMatrixType toDenseMatrix() const 118 { 119 DenseMatrixType res(rows(), cols()); 120 evalToLazy(res); 121 return res; 122 } 123 124 protected: 125 126 void check_coordinates(Index row, Index col) const 127 { 128 EIGEN_ONLY_USED_FOR_DEBUG(row); 129 EIGEN_ONLY_USED_FOR_DEBUG(col); 130 eigen_assert(col>=0 && col<cols() && row>=0 && row<rows()); 131 const int mode = int(Mode) & ~SelfAdjoint; 132 EIGEN_ONLY_USED_FOR_DEBUG(mode); 133 eigen_assert((mode==Upper && col>=row) 134 || (mode==Lower && col<=row) 135 || ((mode==StrictlyUpper || mode==UnitUpper) && col>row) 136 || ((mode==StrictlyLower || mode==UnitLower) && col<row)); 137 } 138 139 #ifdef EIGEN_INTERNAL_DEBUGGING 140 void check_coordinates_internal(Index row, Index col) const 141 { 142 check_coordinates(row, col); 143 } 144 #else 145 void check_coordinates_internal(Index , Index ) const {} 146 #endif 147 148 }; 149 150 /** \class TriangularView 151 * \ingroup Core_Module 152 * 153 * \brief Expression of a triangular part in a matrix 154 * 155 * \param MatrixType the type of the object in which we are taking the triangular part 156 * \param Mode the kind of triangular matrix expression to construct. Can be #Upper, 157 * #Lower, #UnitUpper, #UnitLower, #StrictlyUpper, or #StrictlyLower. 158 * This is in fact a bit field; it must have either #Upper or #Lower, 159 * and additionally it may have #UnitDiag or #ZeroDiag or neither. 160 * 161 * This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular 162 * matrices one should speak of "trapezoid" parts. This class is the return type 163 * of MatrixBase::triangularView() and SparseMatrixBase::triangularView(), and most of the time this is the only way it is used. 164 * 165 * \sa MatrixBase::triangularView() 166 */ 167 namespace internal { 168 template<typename MatrixType, unsigned int _Mode> 169 struct traits<TriangularView<MatrixType, _Mode> > : traits<MatrixType> 170 { 171 typedef typename ref_selector<MatrixType>::non_const_type MatrixTypeNested; 172 typedef typename remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef; 173 typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned; 174 typedef typename MatrixType::PlainObject FullMatrixType; 175 typedef MatrixType ExpressionType; 176 enum { 177 Mode = _Mode, 178 FlagsLvalueBit = is_lvalue<MatrixType>::value ? LvalueBit : 0, 179 Flags = (MatrixTypeNestedCleaned::Flags & (HereditaryBits | FlagsLvalueBit) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) 180 }; 181 }; 182 } 183 184 template<typename _MatrixType, unsigned int _Mode, typename StorageKind> class TriangularViewImpl; 185 186 template<typename _MatrixType, unsigned int _Mode> class TriangularView 187 : public TriangularViewImpl<_MatrixType, _Mode, typename internal::traits<_MatrixType>::StorageKind > 188 { 189 public: 190 191 typedef TriangularViewImpl<_MatrixType, _Mode, typename internal::traits<_MatrixType>::StorageKind > Base; 192 typedef typename internal::traits<TriangularView>::Scalar Scalar; 193 typedef _MatrixType MatrixType; 194 195 protected: 196 typedef typename internal::traits<TriangularView>::MatrixTypeNested MatrixTypeNested; 197 typedef typename internal::traits<TriangularView>::MatrixTypeNestedNonRef MatrixTypeNestedNonRef; 198 199 typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType; 200 201 public: 202 203 typedef typename internal::traits<TriangularView>::StorageKind StorageKind; 204 typedef typename internal::traits<TriangularView>::MatrixTypeNestedCleaned NestedExpression; 205 206 enum { 207 Mode = _Mode, 208 Flags = internal::traits<TriangularView>::Flags, 209 TransposeMode = (Mode & Upper ? Lower : 0) 210 | (Mode & Lower ? Upper : 0) 211 | (Mode & (UnitDiag)) 212 | (Mode & (ZeroDiag)), 213 IsVectorAtCompileTime = false 214 }; 215 216 EIGEN_DEVICE_FUNC 217 explicit inline TriangularView(MatrixType& matrix) : m_matrix(matrix) 218 {} 219 220 using Base::operator=; 221 TriangularView& operator=(const TriangularView &other) 222 { return Base::operator=(other); } 223 224 /** \copydoc EigenBase::rows() */ 225 EIGEN_DEVICE_FUNC 226 inline Index rows() const { return m_matrix.rows(); } 227 /** \copydoc EigenBase::cols() */ 228 EIGEN_DEVICE_FUNC 229 inline Index cols() const { return m_matrix.cols(); } 230 231 /** \returns a const reference to the nested expression */ 232 EIGEN_DEVICE_FUNC 233 const NestedExpression& nestedExpression() const { return m_matrix; } 234 235 /** \returns a reference to the nested expression */ 236 EIGEN_DEVICE_FUNC 237 NestedExpression& nestedExpression() { return m_matrix; } 238 239 typedef TriangularView<const MatrixConjugateReturnType,Mode> ConjugateReturnType; 240 /** \sa MatrixBase::conjugate() const */ 241 EIGEN_DEVICE_FUNC 242 inline const ConjugateReturnType conjugate() const 243 { return ConjugateReturnType(m_matrix.conjugate()); } 244 245 typedef TriangularView<const typename MatrixType::AdjointReturnType,TransposeMode> AdjointReturnType; 246 /** \sa MatrixBase::adjoint() const */ 247 EIGEN_DEVICE_FUNC 248 inline const AdjointReturnType adjoint() const 249 { return AdjointReturnType(m_matrix.adjoint()); } 250 251 typedef TriangularView<typename MatrixType::TransposeReturnType,TransposeMode> TransposeReturnType; 252 /** \sa MatrixBase::transpose() */ 253 EIGEN_DEVICE_FUNC 254 inline TransposeReturnType transpose() 255 { 256 EIGEN_STATIC_ASSERT_LVALUE(MatrixType) 257 typename MatrixType::TransposeReturnType tmp(m_matrix); 258 return TransposeReturnType(tmp); 259 } 260 261 typedef TriangularView<const typename MatrixType::ConstTransposeReturnType,TransposeMode> ConstTransposeReturnType; 262 /** \sa MatrixBase::transpose() const */ 263 EIGEN_DEVICE_FUNC 264 inline const ConstTransposeReturnType transpose() const 265 { 266 return ConstTransposeReturnType(m_matrix.transpose()); 267 } 268 269 template<typename Other> 270 EIGEN_DEVICE_FUNC 271 inline const Solve<TriangularView, Other> 272 solve(const MatrixBase<Other>& other) const 273 { return Solve<TriangularView, Other>(*this, other.derived()); } 274 275 // workaround MSVC ICE 276 #if EIGEN_COMP_MSVC 277 template<int Side, typename Other> 278 EIGEN_DEVICE_FUNC 279 inline const internal::triangular_solve_retval<Side,TriangularView, Other> 280 solve(const MatrixBase<Other>& other) const 281 { return Base::template solve<Side>(other); } 282 #else 283 using Base::solve; 284 #endif 285 286 /** \returns a selfadjoint view of the referenced triangular part which must be either \c #Upper or \c #Lower. 287 * 288 * This is a shortcut for \code this->nestedExpression().selfadjointView<(*this)::Mode>() \endcode 289 * \sa MatrixBase::selfadjointView() */ 290 EIGEN_DEVICE_FUNC 291 SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() 292 { 293 EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR); 294 return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix); 295 } 296 297 /** This is the const version of selfadjointView() */ 298 EIGEN_DEVICE_FUNC 299 const SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() const 300 { 301 EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR); 302 return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix); 303 } 304 305 306 /** \returns the determinant of the triangular matrix 307 * \sa MatrixBase::determinant() */ 308 EIGEN_DEVICE_FUNC 309 Scalar determinant() const 310 { 311 if (Mode & UnitDiag) 312 return 1; 313 else if (Mode & ZeroDiag) 314 return 0; 315 else 316 return m_matrix.diagonal().prod(); 317 } 318 319 protected: 320 321 MatrixTypeNested m_matrix; 322 }; 323 324 /** \ingroup Core_Module 325 * 326 * \brief Base class for a triangular part in a \b dense matrix 327 * 328 * This class is an abstract base class of class TriangularView, and objects of type TriangularViewImpl cannot be instantiated. 329 * It extends class TriangularView with additional methods which available for dense expressions only. 330 * 331 * \sa class TriangularView, MatrixBase::triangularView() 332 */ 333 template<typename _MatrixType, unsigned int _Mode> class TriangularViewImpl<_MatrixType,_Mode,Dense> 334 : public TriangularBase<TriangularView<_MatrixType, _Mode> > 335 { 336 public: 337 338 typedef TriangularView<_MatrixType, _Mode> TriangularViewType; 339 typedef TriangularBase<TriangularViewType> Base; 340 typedef typename internal::traits<TriangularViewType>::Scalar Scalar; 341 342 typedef _MatrixType MatrixType; 343 typedef typename MatrixType::PlainObject DenseMatrixType; 344 typedef DenseMatrixType PlainObject; 345 346 public: 347 using Base::evalToLazy; 348 using Base::derived; 349 350 typedef typename internal::traits<TriangularViewType>::StorageKind StorageKind; 351 352 enum { 353 Mode = _Mode, 354 Flags = internal::traits<TriangularViewType>::Flags 355 }; 356 357 /** \returns the outer-stride of the underlying dense matrix 358 * \sa DenseCoeffsBase::outerStride() */ 359 EIGEN_DEVICE_FUNC 360 inline Index outerStride() const { return derived().nestedExpression().outerStride(); } 361 /** \returns the inner-stride of the underlying dense matrix 362 * \sa DenseCoeffsBase::innerStride() */ 363 EIGEN_DEVICE_FUNC 364 inline Index innerStride() const { return derived().nestedExpression().innerStride(); } 365 366 /** \sa MatrixBase::operator+=() */ 367 template<typename Other> 368 EIGEN_DEVICE_FUNC 369 TriangularViewType& operator+=(const DenseBase<Other>& other) { 370 internal::call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op<Scalar,typename Other::Scalar>()); 371 return derived(); 372 } 373 /** \sa MatrixBase::operator-=() */ 374 template<typename Other> 375 EIGEN_DEVICE_FUNC 376 TriangularViewType& operator-=(const DenseBase<Other>& other) { 377 internal::call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op<Scalar,typename Other::Scalar>()); 378 return derived(); 379 } 380 381 /** \sa MatrixBase::operator*=() */ 382 EIGEN_DEVICE_FUNC 383 TriangularViewType& operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() * other; } 384 /** \sa DenseBase::operator/=() */ 385 EIGEN_DEVICE_FUNC 386 TriangularViewType& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() / other; } 387 388 /** \sa MatrixBase::fill() */ 389 EIGEN_DEVICE_FUNC 390 void fill(const Scalar& value) { setConstant(value); } 391 /** \sa MatrixBase::setConstant() */ 392 EIGEN_DEVICE_FUNC 393 TriangularViewType& setConstant(const Scalar& value) 394 { return *this = MatrixType::Constant(derived().rows(), derived().cols(), value); } 395 /** \sa MatrixBase::setZero() */ 396 EIGEN_DEVICE_FUNC 397 TriangularViewType& setZero() { return setConstant(Scalar(0)); } 398 /** \sa MatrixBase::setOnes() */ 399 EIGEN_DEVICE_FUNC 400 TriangularViewType& setOnes() { return setConstant(Scalar(1)); } 401 402 /** \sa MatrixBase::coeff() 403 * \warning the coordinates must fit into the referenced triangular part 404 */ 405 EIGEN_DEVICE_FUNC 406 inline Scalar coeff(Index row, Index col) const 407 { 408 Base::check_coordinates_internal(row, col); 409 return derived().nestedExpression().coeff(row, col); 410 } 411 412 /** \sa MatrixBase::coeffRef() 413 * \warning the coordinates must fit into the referenced triangular part 414 */ 415 EIGEN_DEVICE_FUNC 416 inline Scalar& coeffRef(Index row, Index col) 417 { 418 EIGEN_STATIC_ASSERT_LVALUE(TriangularViewType); 419 Base::check_coordinates_internal(row, col); 420 return derived().nestedExpression().coeffRef(row, col); 421 } 422 423 /** Assigns a triangular matrix to a triangular part of a dense matrix */ 424 template<typename OtherDerived> 425 EIGEN_DEVICE_FUNC 426 TriangularViewType& operator=(const TriangularBase<OtherDerived>& other); 427 428 /** Shortcut for\code *this = other.other.triangularView<(*this)::Mode>() \endcode */ 429 template<typename OtherDerived> 430 EIGEN_DEVICE_FUNC 431 TriangularViewType& operator=(const MatrixBase<OtherDerived>& other); 432 433 #ifndef EIGEN_PARSED_BY_DOXYGEN 434 EIGEN_DEVICE_FUNC 435 TriangularViewType& operator=(const TriangularViewImpl& other) 436 { return *this = other.derived().nestedExpression(); } 437 438 /** \deprecated */ 439 template<typename OtherDerived> 440 EIGEN_DEVICE_FUNC 441 void lazyAssign(const TriangularBase<OtherDerived>& other); 442 443 /** \deprecated */ 444 template<typename OtherDerived> 445 EIGEN_DEVICE_FUNC 446 void lazyAssign(const MatrixBase<OtherDerived>& other); 447 #endif 448 449 /** Efficient triangular matrix times vector/matrix product */ 450 template<typename OtherDerived> 451 EIGEN_DEVICE_FUNC 452 const Product<TriangularViewType,OtherDerived> 453 operator*(const MatrixBase<OtherDerived>& rhs) const 454 { 455 return Product<TriangularViewType,OtherDerived>(derived(), rhs.derived()); 456 } 457 458 /** Efficient vector/matrix times triangular matrix product */ 459 template<typename OtherDerived> friend 460 EIGEN_DEVICE_FUNC 461 const Product<OtherDerived,TriangularViewType> 462 operator*(const MatrixBase<OtherDerived>& lhs, const TriangularViewImpl& rhs) 463 { 464 return Product<OtherDerived,TriangularViewType>(lhs.derived(),rhs.derived()); 465 } 466 467 /** \returns the product of the inverse of \c *this with \a other, \a *this being triangular. 468 * 469 * This function computes the inverse-matrix matrix product inverse(\c *this) * \a other if 470 * \a Side==OnTheLeft (the default), or the right-inverse-multiply \a other * inverse(\c *this) if 471 * \a Side==OnTheRight. 472 * 473 * Note that the template parameter \c Side can be ommitted, in which case \c Side==OnTheLeft 474 * 475 * The matrix \c *this must be triangular and invertible (i.e., all the coefficients of the 476 * diagonal must be non zero). It works as a forward (resp. backward) substitution if \c *this 477 * is an upper (resp. lower) triangular matrix. 478 * 479 * Example: \include Triangular_solve.cpp 480 * Output: \verbinclude Triangular_solve.out 481 * 482 * This function returns an expression of the inverse-multiply and can works in-place if it is assigned 483 * to the same matrix or vector \a other. 484 * 485 * For users coming from BLAS, this function (and more specifically solveInPlace()) offer 486 * all the operations supported by the \c *TRSV and \c *TRSM BLAS routines. 487 * 488 * \sa TriangularView::solveInPlace() 489 */ 490 template<int Side, typename Other> 491 EIGEN_DEVICE_FUNC 492 inline const internal::triangular_solve_retval<Side,TriangularViewType, Other> 493 solve(const MatrixBase<Other>& other) const; 494 495 /** "in-place" version of TriangularView::solve() where the result is written in \a other 496 * 497 * \warning The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here. 498 * This function will const_cast it, so constness isn't honored here. 499 * 500 * Note that the template parameter \c Side can be ommitted, in which case \c Side==OnTheLeft 501 * 502 * See TriangularView:solve() for the details. 503 */ 504 template<int Side, typename OtherDerived> 505 EIGEN_DEVICE_FUNC 506 void solveInPlace(const MatrixBase<OtherDerived>& other) const; 507 508 template<typename OtherDerived> 509 EIGEN_DEVICE_FUNC 510 void solveInPlace(const MatrixBase<OtherDerived>& other) const 511 { return solveInPlace<OnTheLeft>(other); } 512 513 /** Swaps the coefficients of the common triangular parts of two matrices */ 514 template<typename OtherDerived> 515 EIGEN_DEVICE_FUNC 516 #ifdef EIGEN_PARSED_BY_DOXYGEN 517 void swap(TriangularBase<OtherDerived> &other) 518 #else 519 void swap(TriangularBase<OtherDerived> const & other) 520 #endif 521 { 522 EIGEN_STATIC_ASSERT_LVALUE(OtherDerived); 523 call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>()); 524 } 525 526 /** \deprecated 527 * Shortcut for \code (*this).swap(other.triangularView<(*this)::Mode>()) \endcode */ 528 template<typename OtherDerived> 529 EIGEN_DEVICE_FUNC 530 void swap(MatrixBase<OtherDerived> const & other) 531 { 532 EIGEN_STATIC_ASSERT_LVALUE(OtherDerived); 533 call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>()); 534 } 535 536 template<typename RhsType, typename DstType> 537 EIGEN_DEVICE_FUNC 538 EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const { 539 if(!internal::is_same_dense(dst,rhs)) 540 dst = rhs; 541 this->solveInPlace(dst); 542 } 543 544 template<typename ProductType> 545 EIGEN_DEVICE_FUNC 546 EIGEN_STRONG_INLINE TriangularViewType& _assignProduct(const ProductType& prod, const Scalar& alpha, bool beta); 547 }; 548 549 /*************************************************************************** 550 * Implementation of triangular evaluation/assignment 551 ***************************************************************************/ 552 553 #ifndef EIGEN_PARSED_BY_DOXYGEN 554 // FIXME should we keep that possibility 555 template<typename MatrixType, unsigned int Mode> 556 template<typename OtherDerived> 557 inline TriangularView<MatrixType, Mode>& 558 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const MatrixBase<OtherDerived>& other) 559 { 560 internal::call_assignment_no_alias(derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>()); 561 return derived(); 562 } 563 564 // FIXME should we keep that possibility 565 template<typename MatrixType, unsigned int Mode> 566 template<typename OtherDerived> 567 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const MatrixBase<OtherDerived>& other) 568 { 569 internal::call_assignment_no_alias(derived(), other.template triangularView<Mode>()); 570 } 571 572 573 574 template<typename MatrixType, unsigned int Mode> 575 template<typename OtherDerived> 576 inline TriangularView<MatrixType, Mode>& 577 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const TriangularBase<OtherDerived>& other) 578 { 579 eigen_assert(Mode == int(OtherDerived::Mode)); 580 internal::call_assignment(derived(), other.derived()); 581 return derived(); 582 } 583 584 template<typename MatrixType, unsigned int Mode> 585 template<typename OtherDerived> 586 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const TriangularBase<OtherDerived>& other) 587 { 588 eigen_assert(Mode == int(OtherDerived::Mode)); 589 internal::call_assignment_no_alias(derived(), other.derived()); 590 } 591 #endif 592 593 /*************************************************************************** 594 * Implementation of TriangularBase methods 595 ***************************************************************************/ 596 597 /** Assigns a triangular or selfadjoint matrix to a dense matrix. 598 * If the matrix is triangular, the opposite part is set to zero. */ 599 template<typename Derived> 600 template<typename DenseDerived> 601 void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const 602 { 603 evalToLazy(other.derived()); 604 } 605 606 /*************************************************************************** 607 * Implementation of TriangularView methods 608 ***************************************************************************/ 609 610 /*************************************************************************** 611 * Implementation of MatrixBase methods 612 ***************************************************************************/ 613 614 /** 615 * \returns an expression of a triangular view extracted from the current matrix 616 * 617 * The parameter \a Mode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper, 618 * \c #Lower, \c #StrictlyLower, \c #UnitLower. 619 * 620 * Example: \include MatrixBase_triangularView.cpp 621 * Output: \verbinclude MatrixBase_triangularView.out 622 * 623 * \sa class TriangularView 624 */ 625 template<typename Derived> 626 template<unsigned int Mode> 627 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type 628 MatrixBase<Derived>::triangularView() 629 { 630 return typename TriangularViewReturnType<Mode>::Type(derived()); 631 } 632 633 /** This is the const version of MatrixBase::triangularView() */ 634 template<typename Derived> 635 template<unsigned int Mode> 636 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type 637 MatrixBase<Derived>::triangularView() const 638 { 639 return typename ConstTriangularViewReturnType<Mode>::Type(derived()); 640 } 641 642 /** \returns true if *this is approximately equal to an upper triangular matrix, 643 * within the precision given by \a prec. 644 * 645 * \sa isLowerTriangular() 646 */ 647 template<typename Derived> 648 bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const 649 { 650 RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1); 651 for(Index j = 0; j < cols(); ++j) 652 { 653 Index maxi = numext::mini(j, rows()-1); 654 for(Index i = 0; i <= maxi; ++i) 655 { 656 RealScalar absValue = numext::abs(coeff(i,j)); 657 if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue; 658 } 659 } 660 RealScalar threshold = maxAbsOnUpperPart * prec; 661 for(Index j = 0; j < cols(); ++j) 662 for(Index i = j+1; i < rows(); ++i) 663 if(numext::abs(coeff(i, j)) > threshold) return false; 664 return true; 665 } 666 667 /** \returns true if *this is approximately equal to a lower triangular matrix, 668 * within the precision given by \a prec. 669 * 670 * \sa isUpperTriangular() 671 */ 672 template<typename Derived> 673 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const 674 { 675 RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1); 676 for(Index j = 0; j < cols(); ++j) 677 for(Index i = j; i < rows(); ++i) 678 { 679 RealScalar absValue = numext::abs(coeff(i,j)); 680 if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue; 681 } 682 RealScalar threshold = maxAbsOnLowerPart * prec; 683 for(Index j = 1; j < cols(); ++j) 684 { 685 Index maxi = numext::mini(j, rows()-1); 686 for(Index i = 0; i < maxi; ++i) 687 if(numext::abs(coeff(i, j)) > threshold) return false; 688 } 689 return true; 690 } 691 692 693 /*************************************************************************** 694 **************************************************************************** 695 * Evaluators and Assignment of triangular expressions 696 *************************************************************************** 697 ***************************************************************************/ 698 699 namespace internal { 700 701 702 // TODO currently a triangular expression has the form TriangularView<.,.> 703 // in the future triangular-ness should be defined by the expression traits 704 // such that Transpose<TriangularView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work) 705 template<typename MatrixType, unsigned int Mode> 706 struct evaluator_traits<TriangularView<MatrixType,Mode> > 707 { 708 typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind; 709 typedef typename glue_shapes<typename evaluator_traits<MatrixType>::Shape, TriangularShape>::type Shape; 710 }; 711 712 template<typename MatrixType, unsigned int Mode> 713 struct unary_evaluator<TriangularView<MatrixType,Mode>, IndexBased> 714 : evaluator<typename internal::remove_all<MatrixType>::type> 715 { 716 typedef TriangularView<MatrixType,Mode> XprType; 717 typedef evaluator<typename internal::remove_all<MatrixType>::type> Base; 718 unary_evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {} 719 }; 720 721 // Additional assignment kinds: 722 struct Triangular2Triangular {}; 723 struct Triangular2Dense {}; 724 struct Dense2Triangular {}; 725 726 727 template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop; 728 729 730 /** \internal Specialization of the dense assignment kernel for triangular matrices. 731 * The main difference is that the triangular, diagonal, and opposite parts are processed through three different functions. 732 * \tparam UpLo must be either Lower or Upper 733 * \tparam Mode must be either 0, UnitDiag, ZeroDiag, or SelfAdjoint 734 */ 735 template<int UpLo, int Mode, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version = Specialized> 736 class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> 737 { 738 protected: 739 typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base; 740 typedef typename Base::DstXprType DstXprType; 741 typedef typename Base::SrcXprType SrcXprType; 742 using Base::m_dst; 743 using Base::m_src; 744 using Base::m_functor; 745 public: 746 747 typedef typename Base::DstEvaluatorType DstEvaluatorType; 748 typedef typename Base::SrcEvaluatorType SrcEvaluatorType; 749 typedef typename Base::Scalar Scalar; 750 typedef typename Base::AssignmentTraits AssignmentTraits; 751 752 753 EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr) 754 : Base(dst, src, func, dstExpr) 755 {} 756 757 #ifdef EIGEN_INTERNAL_DEBUGGING 758 EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col) 759 { 760 eigen_internal_assert(row!=col); 761 Base::assignCoeff(row,col); 762 } 763 #else 764 using Base::assignCoeff; 765 #endif 766 767 EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id) 768 { 769 if(Mode==UnitDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(1)); 770 else if(Mode==ZeroDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(0)); 771 else if(Mode==0) Base::assignCoeff(id,id); 772 } 773 774 EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index row, Index col) 775 { 776 eigen_internal_assert(row!=col); 777 if(SetOpposite) 778 m_functor.assignCoeff(m_dst.coeffRef(row,col), Scalar(0)); 779 } 780 }; 781 782 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType, typename Functor> 783 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE 784 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src, const Functor &func) 785 { 786 typedef evaluator<DstXprType> DstEvaluatorType; 787 typedef evaluator<SrcXprType> SrcEvaluatorType; 788 789 SrcEvaluatorType srcEvaluator(src); 790 791 Index dstRows = src.rows(); 792 Index dstCols = src.cols(); 793 if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) 794 dst.resize(dstRows, dstCols); 795 DstEvaluatorType dstEvaluator(dst); 796 797 typedef triangular_dense_assignment_kernel< Mode&(Lower|Upper),Mode&(UnitDiag|ZeroDiag|SelfAdjoint),SetOpposite, 798 DstEvaluatorType,SrcEvaluatorType,Functor> Kernel; 799 Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived()); 800 801 enum { 802 unroll = DstXprType::SizeAtCompileTime != Dynamic 803 && SrcEvaluatorType::CoeffReadCost < HugeCost 804 && DstXprType::SizeAtCompileTime * (DstEvaluatorType::CoeffReadCost+SrcEvaluatorType::CoeffReadCost) / 2 <= EIGEN_UNROLLING_LIMIT 805 }; 806 807 triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, SetOpposite>::run(kernel); 808 } 809 810 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType> 811 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE 812 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src) 813 { 814 call_triangular_assignment_loop<Mode,SetOpposite>(dst, src, internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>()); 815 } 816 817 template<> struct AssignmentKind<TriangularShape,TriangularShape> { typedef Triangular2Triangular Kind; }; 818 template<> struct AssignmentKind<DenseShape,TriangularShape> { typedef Triangular2Dense Kind; }; 819 template<> struct AssignmentKind<TriangularShape,DenseShape> { typedef Dense2Triangular Kind; }; 820 821 822 template< typename DstXprType, typename SrcXprType, typename Functor> 823 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Triangular> 824 { 825 EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) 826 { 827 eigen_assert(int(DstXprType::Mode) == int(SrcXprType::Mode)); 828 829 call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func); 830 } 831 }; 832 833 template< typename DstXprType, typename SrcXprType, typename Functor> 834 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Dense> 835 { 836 EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) 837 { 838 call_triangular_assignment_loop<SrcXprType::Mode, (SrcXprType::Mode&SelfAdjoint)==0>(dst, src, func); 839 } 840 }; 841 842 template< typename DstXprType, typename SrcXprType, typename Functor> 843 struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular> 844 { 845 EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) 846 { 847 call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func); 848 } 849 }; 850 851 852 template<typename Kernel, unsigned int Mode, int UnrollCount, bool SetOpposite> 853 struct triangular_assignment_loop 854 { 855 // FIXME: this is not very clean, perhaps this information should be provided by the kernel? 856 typedef typename Kernel::DstEvaluatorType DstEvaluatorType; 857 typedef typename DstEvaluatorType::XprType DstXprType; 858 859 enum { 860 col = (UnrollCount-1) / DstXprType::RowsAtCompileTime, 861 row = (UnrollCount-1) % DstXprType::RowsAtCompileTime 862 }; 863 864 typedef typename Kernel::Scalar Scalar; 865 866 EIGEN_DEVICE_FUNC 867 static inline void run(Kernel &kernel) 868 { 869 triangular_assignment_loop<Kernel, Mode, UnrollCount-1, SetOpposite>::run(kernel); 870 871 if(row==col) 872 kernel.assignDiagonalCoeff(row); 873 else if( ((Mode&Lower) && row>col) || ((Mode&Upper) && row<col) ) 874 kernel.assignCoeff(row,col); 875 else if(SetOpposite) 876 kernel.assignOppositeCoeff(row,col); 877 } 878 }; 879 880 // prevent buggy user code from causing an infinite recursion 881 template<typename Kernel, unsigned int Mode, bool SetOpposite> 882 struct triangular_assignment_loop<Kernel, Mode, 0, SetOpposite> 883 { 884 EIGEN_DEVICE_FUNC 885 static inline void run(Kernel &) {} 886 }; 887 888 889 890 // TODO: experiment with a recursive assignment procedure splitting the current 891 // triangular part into one rectangular and two triangular parts. 892 893 894 template<typename Kernel, unsigned int Mode, bool SetOpposite> 895 struct triangular_assignment_loop<Kernel, Mode, Dynamic, SetOpposite> 896 { 897 typedef typename Kernel::Scalar Scalar; 898 EIGEN_DEVICE_FUNC 899 static inline void run(Kernel &kernel) 900 { 901 for(Index j = 0; j < kernel.cols(); ++j) 902 { 903 Index maxi = numext::mini(j, kernel.rows()); 904 Index i = 0; 905 if (((Mode&Lower) && SetOpposite) || (Mode&Upper)) 906 { 907 for(; i < maxi; ++i) 908 if(Mode&Upper) kernel.assignCoeff(i, j); 909 else kernel.assignOppositeCoeff(i, j); 910 } 911 else 912 i = maxi; 913 914 if(i<kernel.rows()) // then i==j 915 kernel.assignDiagonalCoeff(i++); 916 917 if (((Mode&Upper) && SetOpposite) || (Mode&Lower)) 918 { 919 for(; i < kernel.rows(); ++i) 920 if(Mode&Lower) kernel.assignCoeff(i, j); 921 else kernel.assignOppositeCoeff(i, j); 922 } 923 } 924 } 925 }; 926 927 } // end namespace internal 928 929 /** Assigns a triangular or selfadjoint matrix to a dense matrix. 930 * If the matrix is triangular, the opposite part is set to zero. */ 931 template<typename Derived> 932 template<typename DenseDerived> 933 void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const 934 { 935 other.derived().resize(this->rows(), this->cols()); 936 internal::call_triangular_assignment_loop<Derived::Mode,(Derived::Mode&SelfAdjoint)==0 /* SetOpposite */>(other.derived(), derived().nestedExpression()); 937 } 938 939 namespace internal { 940 941 // Triangular = Product 942 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar> 943 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular> 944 { 945 typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType; 946 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,typename SrcXprType::Scalar> &) 947 { 948 Index dstRows = src.rows(); 949 Index dstCols = src.cols(); 950 if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) 951 dst.resize(dstRows, dstCols); 952 953 dst._assignProduct(src, 1, 0); 954 } 955 }; 956 957 // Triangular += Product 958 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar> 959 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::add_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular> 960 { 961 typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType; 962 static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar,typename SrcXprType::Scalar> &) 963 { 964 dst._assignProduct(src, 1, 1); 965 } 966 }; 967 968 // Triangular -= Product 969 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar> 970 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::sub_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular> 971 { 972 typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType; 973 static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,typename SrcXprType::Scalar> &) 974 { 975 dst._assignProduct(src, -1, 1); 976 } 977 }; 978 979 } // end namespace internal 980 981 } // end namespace Eigen 982 983 #endif // EIGEN_TRIANGULARMATRIX_H 984
