00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef EIGEN_TRIANGULARMATRIX_H
00027 #define EIGEN_TRIANGULARMATRIX_H
00028
00029 namespace Eigen {
00030
00031 namespace internal {
00032
00033 template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval;
00034
00035 }
00036
00044 template<typename Derived> class TriangularBase : public EigenBase<Derived>
00045 {
00046 public:
00047
00048 enum {
00049 Mode = internal::traits<Derived>::Mode,
00050 CoeffReadCost = internal::traits<Derived>::CoeffReadCost,
00051 RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
00052 ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
00053 MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
00054 MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime
00055 };
00056 typedef typename internal::traits<Derived>::Scalar Scalar;
00057 typedef typename internal::traits<Derived>::StorageKind StorageKind;
00058 typedef typename internal::traits<Derived>::Index Index;
00059 typedef typename internal::traits<Derived>::DenseMatrixType DenseMatrixType;
00060 typedef DenseMatrixType DenseType;
00061
00062 inline TriangularBase() { eigen_assert(!((Mode&UnitDiag) && (Mode&ZeroDiag))); }
00063
00064 inline Index rows() const { return derived().rows(); }
00065 inline Index cols() const { return derived().cols(); }
00066 inline Index outerStride() const { return derived().outerStride(); }
00067 inline Index innerStride() const { return derived().innerStride(); }
00068
00069 inline Scalar coeff(Index row, Index col) const { return derived().coeff(row,col); }
00070 inline Scalar& coeffRef(Index row, Index col) { return derived().coeffRef(row,col); }
00071
00074 template<typename Other>
00075 EIGEN_STRONG_INLINE void copyCoeff(Index row, Index col, Other& other)
00076 {
00077 derived().coeffRef(row, col) = other.coeff(row, col);
00078 }
00079
00080 inline Scalar operator()(Index row, Index col) const
00081 {
00082 check_coordinates(row, col);
00083 return coeff(row,col);
00084 }
00085 inline Scalar& operator()(Index row, Index col)
00086 {
00087 check_coordinates(row, col);
00088 return coeffRef(row,col);
00089 }
00090
00091 #ifndef EIGEN_PARSED_BY_DOXYGEN
00092 inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
00093 inline Derived& derived() { return *static_cast<Derived*>(this); }
00094 #endif // not EIGEN_PARSED_BY_DOXYGEN
00095
00096 template<typename DenseDerived>
00097 void evalTo(MatrixBase<DenseDerived> &other) const;
00098 template<typename DenseDerived>
00099 void evalToLazy(MatrixBase<DenseDerived> &other) const;
00100
00101 DenseMatrixType toDenseMatrix() const
00102 {
00103 DenseMatrixType res(rows(), cols());
00104 evalToLazy(res);
00105 return res;
00106 }
00107
00108 protected:
00109
00110 void check_coordinates(Index row, Index col) const
00111 {
00112 EIGEN_ONLY_USED_FOR_DEBUG(row);
00113 EIGEN_ONLY_USED_FOR_DEBUG(col);
00114 eigen_assert(col>=0 && col<cols() && row>=0 && row<rows());
00115 const int mode = int(Mode) & ~SelfAdjoint;
00116 EIGEN_ONLY_USED_FOR_DEBUG(mode);
00117 eigen_assert((mode==Upper && col>=row)
00118 || (mode==Lower && col<=row)
00119 || ((mode==StrictlyUpper || mode==UnitUpper) && col>row)
00120 || ((mode==StrictlyLower || mode==UnitLower) && col<row));
00121 }
00122
00123 #ifdef EIGEN_INTERNAL_DEBUGGING
00124 void check_coordinates_internal(Index row, Index col) const
00125 {
00126 check_coordinates(row, col);
00127 }
00128 #else
00129 void check_coordinates_internal(Index , Index ) const {}
00130 #endif
00131
00132 };
00133
00151 namespace internal {
00152 template<typename MatrixType, unsigned int _Mode>
00153 struct traits<TriangularView<MatrixType, _Mode> > : traits<MatrixType>
00154 {
00155 typedef typename nested<MatrixType>::type MatrixTypeNested;
00156 typedef typename remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef;
00157 typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
00158 typedef MatrixType ExpressionType;
00159 typedef typename MatrixType::PlainObject DenseMatrixType;
00160 enum {
00161 Mode = _Mode,
00162 Flags = (MatrixTypeNestedCleaned::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) | Mode,
00163 CoeffReadCost = MatrixTypeNestedCleaned::CoeffReadCost
00164 };
00165 };
00166 }
00167
00168 template<int Mode, bool LhsIsTriangular,
00169 typename Lhs, bool LhsIsVector,
00170 typename Rhs, bool RhsIsVector>
00171 struct TriangularProduct;
00172
00173 template<typename _MatrixType, unsigned int _Mode> class TriangularView
00174 : public TriangularBase<TriangularView<_MatrixType, _Mode> >
00175 {
00176 public:
00177
00178 typedef TriangularBase<TriangularView> Base;
00179 typedef typename internal::traits<TriangularView>::Scalar Scalar;
00180
00181 typedef _MatrixType MatrixType;
00182 typedef typename internal::traits<TriangularView>::DenseMatrixType DenseMatrixType;
00183 typedef DenseMatrixType PlainObject;
00184
00185 protected:
00186 typedef typename internal::traits<TriangularView>::MatrixTypeNested MatrixTypeNested;
00187 typedef typename internal::traits<TriangularView>::MatrixTypeNestedNonRef MatrixTypeNestedNonRef;
00188 typedef typename internal::traits<TriangularView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned;
00189
00190 typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType;
00191
00192 public:
00193 using Base::evalToLazy;
00194
00195
00196 typedef typename internal::traits<TriangularView>::StorageKind StorageKind;
00197 typedef typename internal::traits<TriangularView>::Index Index;
00198
00199 enum {
00200 Mode = _Mode,
00201 TransposeMode = (Mode & Upper ? Lower : 0)
00202 | (Mode & Lower ? Upper : 0)
00203 | (Mode & (UnitDiag))
00204 | (Mode & (ZeroDiag))
00205 };
00206
00207 inline TriangularView(const MatrixType& matrix) : m_matrix(matrix)
00208 {}
00209
00210 inline Index rows() const { return m_matrix.rows(); }
00211 inline Index cols() const { return m_matrix.cols(); }
00212 inline Index outerStride() const { return m_matrix.outerStride(); }
00213 inline Index innerStride() const { return m_matrix.innerStride(); }
00214
00216 template<typename Other> TriangularView& operator+=(const DenseBase<Other>& other) { return *this = m_matrix + other.derived(); }
00218 template<typename Other> TriangularView& operator-=(const DenseBase<Other>& other) { return *this = m_matrix - other.derived(); }
00220 TriangularView& operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix * other; }
00222 TriangularView& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix / other; }
00223
00225 void fill(const Scalar& value) { setConstant(value); }
00227 TriangularView& setConstant(const Scalar& value)
00228 { return *this = MatrixType::Constant(rows(), cols(), value); }
00230 TriangularView& setZero() { return setConstant(Scalar(0)); }
00232 TriangularView& setOnes() { return setConstant(Scalar(1)); }
00233
00237 inline Scalar coeff(Index row, Index col) const
00238 {
00239 Base::check_coordinates_internal(row, col);
00240 return m_matrix.coeff(row, col);
00241 }
00242
00246 inline Scalar& coeffRef(Index row, Index col)
00247 {
00248 Base::check_coordinates_internal(row, col);
00249 return m_matrix.const_cast_derived().coeffRef(row, col);
00250 }
00251
00252 const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; }
00253 MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); }
00254
00256 template<typename OtherDerived>
00257 TriangularView& operator=(const TriangularBase<OtherDerived>& other);
00258
00259 template<typename OtherDerived>
00260 TriangularView& operator=(const MatrixBase<OtherDerived>& other);
00261
00262 TriangularView& operator=(const TriangularView& other)
00263 { return *this = other.nestedExpression(); }
00264
00265 template<typename OtherDerived>
00266 void lazyAssign(const TriangularBase<OtherDerived>& other);
00267
00268 template<typename OtherDerived>
00269 void lazyAssign(const MatrixBase<OtherDerived>& other);
00270
00272 inline TriangularView<MatrixConjugateReturnType,Mode> conjugate()
00273 { return m_matrix.conjugate(); }
00275 inline const TriangularView<MatrixConjugateReturnType,Mode> conjugate() const
00276 { return m_matrix.conjugate(); }
00277
00279 inline const TriangularView<const typename MatrixType::AdjointReturnType,TransposeMode> adjoint() const
00280 { return m_matrix.adjoint(); }
00281
00283 inline TriangularView<Transpose<MatrixType>,TransposeMode> transpose()
00284 {
00285 EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
00286 return m_matrix.const_cast_derived().transpose();
00287 }
00289 inline const TriangularView<Transpose<MatrixType>,TransposeMode> transpose() const
00290 {
00291 return m_matrix.transpose();
00292 }
00293
00295 template<typename OtherDerived>
00296 TriangularProduct<Mode,true,MatrixType,false,OtherDerived, OtherDerived::IsVectorAtCompileTime>
00297 operator*(const MatrixBase<OtherDerived>& rhs) const
00298 {
00299 return TriangularProduct
00300 <Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime>
00301 (m_matrix, rhs.derived());
00302 }
00303
00305 template<typename OtherDerived> friend
00306 TriangularProduct<Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false>
00307 operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs)
00308 {
00309 return TriangularProduct
00310 <Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false>
00311 (lhs.derived(),rhs.m_matrix);
00312 }
00313
00314 #ifdef EIGEN2_SUPPORT
00315 template<typename OtherDerived>
00316 struct eigen2_product_return_type
00317 {
00318 typedef typename TriangularView<MatrixType,Mode>::DenseMatrixType DenseMatrixType;
00319 typedef typename OtherDerived::PlainObject::DenseType OtherPlainObject;
00320 typedef typename ProductReturnType<DenseMatrixType, OtherPlainObject>::Type ProdRetType;
00321 typedef typename ProdRetType::PlainObject type;
00322 };
00323 template<typename OtherDerived>
00324 const typename eigen2_product_return_type<OtherDerived>::type
00325 operator*(const EigenBase<OtherDerived>& rhs) const
00326 {
00327 typename OtherDerived::PlainObject::DenseType rhsPlainObject;
00328 rhs.evalTo(rhsPlainObject);
00329 return this->toDenseMatrix() * rhsPlainObject;
00330 }
00331 template<typename OtherMatrixType>
00332 bool isApprox(const TriangularView<OtherMatrixType, Mode>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const
00333 {
00334 return this->toDenseMatrix().isApprox(other.toDenseMatrix(), precision);
00335 }
00336 template<typename OtherDerived>
00337 bool isApprox(const MatrixBase<OtherDerived>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const
00338 {
00339 return this->toDenseMatrix().isApprox(other, precision);
00340 }
00341 #endif // EIGEN2_SUPPORT
00342
00343 template<int Side, typename Other>
00344 inline const internal::triangular_solve_retval<Side,TriangularView, Other>
00345 solve(const MatrixBase<Other>& other) const;
00346
00347 template<int Side, typename OtherDerived>
00348 void solveInPlace(const MatrixBase<OtherDerived>& other) const;
00349
00350 template<typename Other>
00351 inline const internal::triangular_solve_retval<OnTheLeft,TriangularView, Other>
00352 solve(const MatrixBase<Other>& other) const
00353 { return solve<OnTheLeft>(other); }
00354
00355 template<typename OtherDerived>
00356 void solveInPlace(const MatrixBase<OtherDerived>& other) const
00357 { return solveInPlace<OnTheLeft>(other); }
00358
00359 const SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() const
00360 {
00361 EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
00362 return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
00363 }
00364 SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView()
00365 {
00366 EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
00367 return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
00368 }
00369
00370 template<typename OtherDerived>
00371 void swap(TriangularBase<OtherDerived> const & other)
00372 {
00373 TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived());
00374 }
00375
00376 template<typename OtherDerived>
00377 void swap(MatrixBase<OtherDerived> const & other)
00378 {
00379 SwapWrapper<MatrixType> swaper(const_cast<MatrixType&>(m_matrix));
00380 TriangularView<SwapWrapper<MatrixType>,Mode>(swaper).lazyAssign(other.derived());
00381 }
00382
00383 Scalar determinant() const
00384 {
00385 if (Mode & UnitDiag)
00386 return 1;
00387 else if (Mode & ZeroDiag)
00388 return 0;
00389 else
00390 return m_matrix.diagonal().prod();
00391 }
00392
00393
00394 template<typename ProductDerived, typename Lhs, typename Rhs>
00395 EIGEN_STRONG_INLINE TriangularView& operator=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
00396 {
00397 setZero();
00398 return assignProduct(other,1);
00399 }
00400
00401 template<typename ProductDerived, typename Lhs, typename Rhs>
00402 EIGEN_STRONG_INLINE TriangularView& operator+=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
00403 {
00404 return assignProduct(other,1);
00405 }
00406
00407 template<typename ProductDerived, typename Lhs, typename Rhs>
00408 EIGEN_STRONG_INLINE TriangularView& operator-=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
00409 {
00410 return assignProduct(other,-1);
00411 }
00412
00413
00414 template<typename ProductDerived>
00415 EIGEN_STRONG_INLINE TriangularView& operator=(const ScaledProduct<ProductDerived>& other)
00416 {
00417 setZero();
00418 return assignProduct(other,other.alpha());
00419 }
00420
00421 template<typename ProductDerived>
00422 EIGEN_STRONG_INLINE TriangularView& operator+=(const ScaledProduct<ProductDerived>& other)
00423 {
00424 return assignProduct(other,other.alpha());
00425 }
00426
00427 template<typename ProductDerived>
00428 EIGEN_STRONG_INLINE TriangularView& operator-=(const ScaledProduct<ProductDerived>& other)
00429 {
00430 return assignProduct(other,-other.alpha());
00431 }
00432
00433 protected:
00434
00435 template<typename ProductDerived, typename Lhs, typename Rhs>
00436 EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase<ProductDerived, Lhs,Rhs>& prod, const Scalar& alpha);
00437
00438 MatrixTypeNested m_matrix;
00439 };
00440
00441
00442
00443
00444
00445 namespace internal {
00446
00447 template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount, bool ClearOpposite>
00448 struct triangular_assignment_selector
00449 {
00450 enum {
00451 col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
00452 row = (UnrollCount-1) % Derived1::RowsAtCompileTime
00453 };
00454
00455 typedef typename Derived1::Scalar Scalar;
00456
00457 static inline void run(Derived1 &dst, const Derived2 &src)
00458 {
00459 triangular_assignment_selector<Derived1, Derived2, Mode, UnrollCount-1, ClearOpposite>::run(dst, src);
00460
00461 eigen_assert( Mode == Upper || Mode == Lower
00462 || Mode == StrictlyUpper || Mode == StrictlyLower
00463 || Mode == UnitUpper || Mode == UnitLower);
00464 if((Mode == Upper && row <= col)
00465 || (Mode == Lower && row >= col)
00466 || (Mode == StrictlyUpper && row < col)
00467 || (Mode == StrictlyLower && row > col)
00468 || (Mode == UnitUpper && row < col)
00469 || (Mode == UnitLower && row > col))
00470 dst.copyCoeff(row, col, src);
00471 else if(ClearOpposite)
00472 {
00473 if (Mode&UnitDiag && row==col)
00474 dst.coeffRef(row, col) = Scalar(1);
00475 else
00476 dst.coeffRef(row, col) = Scalar(0);
00477 }
00478 }
00479 };
00480
00481
00482 template<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite>
00483 struct triangular_assignment_selector<Derived1, Derived2, Mode, 0, ClearOpposite>
00484 {
00485 static inline void run(Derived1 &, const Derived2 &) {}
00486 };
00487
00488 template<typename Derived1, typename Derived2, bool ClearOpposite>
00489 struct triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearOpposite>
00490 {
00491 typedef typename Derived1::Index Index;
00492 typedef typename Derived1::Scalar Scalar;
00493 static inline void run(Derived1 &dst, const Derived2 &src)
00494 {
00495 for(Index j = 0; j < dst.cols(); ++j)
00496 {
00497 Index maxi = (std::min)(j, dst.rows()-1);
00498 for(Index i = 0; i <= maxi; ++i)
00499 dst.copyCoeff(i, j, src);
00500 if (ClearOpposite)
00501 for(Index i = maxi+1; i < dst.rows(); ++i)
00502 dst.coeffRef(i, j) = Scalar(0);
00503 }
00504 }
00505 };
00506
00507 template<typename Derived1, typename Derived2, bool ClearOpposite>
00508 struct triangular_assignment_selector<Derived1, Derived2, Lower, Dynamic, ClearOpposite>
00509 {
00510 typedef typename Derived1::Index Index;
00511 static inline void run(Derived1 &dst, const Derived2 &src)
00512 {
00513 for(Index j = 0; j < dst.cols(); ++j)
00514 {
00515 for(Index i = j; i < dst.rows(); ++i)
00516 dst.copyCoeff(i, j, src);
00517 Index maxi = (std::min)(j, dst.rows());
00518 if (ClearOpposite)
00519 for(Index i = 0; i < maxi; ++i)
00520 dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
00521 }
00522 }
00523 };
00524
00525 template<typename Derived1, typename Derived2, bool ClearOpposite>
00526 struct triangular_assignment_selector<Derived1, Derived2, StrictlyUpper, Dynamic, ClearOpposite>
00527 {
00528 typedef typename Derived1::Index Index;
00529 static inline void run(Derived1 &dst, const Derived2 &src)
00530 {
00531 for(Index j = 0; j < dst.cols(); ++j)
00532 {
00533 Index maxi = (std::min)(j, dst.rows());
00534 for(Index i = 0; i < maxi; ++i)
00535 dst.copyCoeff(i, j, src);
00536 if (ClearOpposite)
00537 for(Index i = maxi; i < dst.rows(); ++i)
00538 dst.coeffRef(i, j) = 0;
00539 }
00540 }
00541 };
00542
00543 template<typename Derived1, typename Derived2, bool ClearOpposite>
00544 struct triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic, ClearOpposite>
00545 {
00546 typedef typename Derived1::Index Index;
00547 static inline void run(Derived1 &dst, const Derived2 &src)
00548 {
00549 for(Index j = 0; j < dst.cols(); ++j)
00550 {
00551 for(Index i = j+1; i < dst.rows(); ++i)
00552 dst.copyCoeff(i, j, src);
00553 Index maxi = (std::min)(j, dst.rows()-1);
00554 if (ClearOpposite)
00555 for(Index i = 0; i <= maxi; ++i)
00556 dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
00557 }
00558 }
00559 };
00560
00561 template<typename Derived1, typename Derived2, bool ClearOpposite>
00562 struct triangular_assignment_selector<Derived1, Derived2, UnitUpper, Dynamic, ClearOpposite>
00563 {
00564 typedef typename Derived1::Index Index;
00565 static inline void run(Derived1 &dst, const Derived2 &src)
00566 {
00567 for(Index j = 0; j < dst.cols(); ++j)
00568 {
00569 Index maxi = (std::min)(j, dst.rows());
00570 for(Index i = 0; i < maxi; ++i)
00571 dst.copyCoeff(i, j, src);
00572 if (ClearOpposite)
00573 {
00574 for(Index i = maxi+1; i < dst.rows(); ++i)
00575 dst.coeffRef(i, j) = 0;
00576 }
00577 }
00578 dst.diagonal().setOnes();
00579 }
00580 };
00581 template<typename Derived1, typename Derived2, bool ClearOpposite>
00582 struct triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, ClearOpposite>
00583 {
00584 typedef typename Derived1::Index Index;
00585 static inline void run(Derived1 &dst, const Derived2 &src)
00586 {
00587 for(Index j = 0; j < dst.cols(); ++j)
00588 {
00589 Index maxi = (std::min)(j, dst.rows());
00590 for(Index i = maxi+1; i < dst.rows(); ++i)
00591 dst.copyCoeff(i, j, src);
00592 if (ClearOpposite)
00593 {
00594 for(Index i = 0; i < maxi; ++i)
00595 dst.coeffRef(i, j) = 0;
00596 }
00597 }
00598 dst.diagonal().setOnes();
00599 }
00600 };
00601
00602 }
00603
00604
00605 template<typename MatrixType, unsigned int Mode>
00606 template<typename OtherDerived>
00607 inline TriangularView<MatrixType, Mode>&
00608 TriangularView<MatrixType, Mode>::operator=(const MatrixBase<OtherDerived>& other)
00609 {
00610 if(OtherDerived::Flags & EvalBeforeAssigningBit)
00611 {
00612 typename internal::plain_matrix_type<OtherDerived>::type other_evaluated(other.rows(), other.cols());
00613 other_evaluated.template triangularView<Mode>().lazyAssign(other.derived());
00614 lazyAssign(other_evaluated);
00615 }
00616 else
00617 lazyAssign(other.derived());
00618 return *this;
00619 }
00620
00621
00622 template<typename MatrixType, unsigned int Mode>
00623 template<typename OtherDerived>
00624 void TriangularView<MatrixType, Mode>::lazyAssign(const MatrixBase<OtherDerived>& other)
00625 {
00626 enum {
00627 unroll = MatrixType::SizeAtCompileTime != Dynamic
00628 && internal::traits<OtherDerived>::CoeffReadCost != Dynamic
00629 && MatrixType::SizeAtCompileTime*internal::traits<OtherDerived>::CoeffReadCost/2 <= EIGEN_UNROLLING_LIMIT
00630 };
00631 eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
00632
00633 internal::triangular_assignment_selector
00634 <MatrixType, OtherDerived, int(Mode),
00635 unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic,
00636 false
00637 >::run(m_matrix.const_cast_derived(), other.derived());
00638 }
00639
00640
00641
00642 template<typename MatrixType, unsigned int Mode>
00643 template<typename OtherDerived>
00644 inline TriangularView<MatrixType, Mode>&
00645 TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>& other)
00646 {
00647 eigen_assert(Mode == int(OtherDerived::Mode));
00648 if(internal::traits<OtherDerived>::Flags & EvalBeforeAssigningBit)
00649 {
00650 typename OtherDerived::DenseMatrixType other_evaluated(other.rows(), other.cols());
00651 other_evaluated.template triangularView<Mode>().lazyAssign(other.derived().nestedExpression());
00652 lazyAssign(other_evaluated);
00653 }
00654 else
00655 lazyAssign(other.derived().nestedExpression());
00656 return *this;
00657 }
00658
00659 template<typename MatrixType, unsigned int Mode>
00660 template<typename OtherDerived>
00661 void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDerived>& other)
00662 {
00663 enum {
00664 unroll = MatrixType::SizeAtCompileTime != Dynamic
00665 && internal::traits<OtherDerived>::CoeffReadCost != Dynamic
00666 && MatrixType::SizeAtCompileTime * internal::traits<OtherDerived>::CoeffReadCost / 2
00667 <= EIGEN_UNROLLING_LIMIT
00668 };
00669 eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
00670
00671 internal::triangular_assignment_selector
00672 <MatrixType, OtherDerived, int(Mode),
00673 unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic,
00674 false
00675 >::run(m_matrix.const_cast_derived(), other.derived().nestedExpression());
00676 }
00677
00678
00679
00680
00681
00684 template<typename Derived>
00685 template<typename DenseDerived>
00686 void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
00687 {
00688 if(internal::traits<Derived>::Flags & EvalBeforeAssigningBit)
00689 {
00690 typename internal::plain_matrix_type<Derived>::type other_evaluated(rows(), cols());
00691 evalToLazy(other_evaluated);
00692 other.derived().swap(other_evaluated);
00693 }
00694 else
00695 evalToLazy(other.derived());
00696 }
00697
00700 template<typename Derived>
00701 template<typename DenseDerived>
00702 void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
00703 {
00704 enum {
00705 unroll = DenseDerived::SizeAtCompileTime != Dynamic
00706 && internal::traits<Derived>::CoeffReadCost != Dynamic
00707 && DenseDerived::SizeAtCompileTime * internal::traits<Derived>::CoeffReadCost / 2
00708 <= EIGEN_UNROLLING_LIMIT
00709 };
00710 other.derived().resize(this->rows(), this->cols());
00711
00712 internal::triangular_assignment_selector
00713 <DenseDerived, typename internal::traits<Derived>::MatrixTypeNestedCleaned, Derived::Mode,
00714 unroll ? int(DenseDerived::SizeAtCompileTime) : Dynamic,
00715 true
00716 >::run(other.derived(), derived().nestedExpression());
00717 }
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727 #ifdef EIGEN2_SUPPORT
00728
00729
00730
00731 namespace internal {
00732 template<typename MatrixType, unsigned int Mode>
00733 struct eigen2_part_return_type
00734 {
00735 typedef TriangularView<MatrixType, Mode> type;
00736 };
00737
00738 template<typename MatrixType>
00739 struct eigen2_part_return_type<MatrixType, SelfAdjoint>
00740 {
00741 typedef SelfAdjointView<MatrixType, Upper> type;
00742 };
00743 }
00744
00746 template<typename Derived>
00747 template<unsigned int Mode>
00748 const typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part() const
00749 {
00750 return derived();
00751 }
00752
00754 template<typename Derived>
00755 template<unsigned int Mode>
00756 typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part()
00757 {
00758 return derived();
00759 }
00760 #endif
00761
00773 template<typename Derived>
00774 template<unsigned int Mode>
00775 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
00776 MatrixBase<Derived>::triangularView()
00777 {
00778 return derived();
00779 }
00780
00782 template<typename Derived>
00783 template<unsigned int Mode>
00784 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
00785 MatrixBase<Derived>::triangularView() const
00786 {
00787 return derived();
00788 }
00789
00795 template<typename Derived>
00796 bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const
00797 {
00798 RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
00799 for(Index j = 0; j < cols(); ++j)
00800 {
00801 Index maxi = (std::min)(j, rows()-1);
00802 for(Index i = 0; i <= maxi; ++i)
00803 {
00804 RealScalar absValue = internal::abs(coeff(i,j));
00805 if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
00806 }
00807 }
00808 RealScalar threshold = maxAbsOnUpperPart * prec;
00809 for(Index j = 0; j < cols(); ++j)
00810 for(Index i = j+1; i < rows(); ++i)
00811 if(internal::abs(coeff(i, j)) > threshold) return false;
00812 return true;
00813 }
00814
00820 template<typename Derived>
00821 bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const
00822 {
00823 RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
00824 for(Index j = 0; j < cols(); ++j)
00825 for(Index i = j; i < rows(); ++i)
00826 {
00827 RealScalar absValue = internal::abs(coeff(i,j));
00828 if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
00829 }
00830 RealScalar threshold = maxAbsOnLowerPart * prec;
00831 for(Index j = 1; j < cols(); ++j)
00832 {
00833 Index maxi = (std::min)(j, rows()-1);
00834 for(Index i = 0; i < maxi; ++i)
00835 if(internal::abs(coeff(i, j)) > threshold) return false;
00836 }
00837 return true;
00838 }
00839
00840 }
00841
00842 #endif // EIGEN_TRIANGULARMATRIX_H