• Home
  • History
  • Annotate
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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::ColsAtCompileTime==1>
282     operator*(const MatrixBase<OtherDerived>& rhs) const
283     {
284       return TriangularProduct
285               <Mode, true, MatrixType, false, OtherDerived, OtherDerived::ColsAtCompileTime==1>
286               (m_matrix, rhs.derived());
287     }
288 
289     /** Efficient vector/matrix times triangular matrix product */
290     template<typename OtherDerived> friend
291     TriangularProduct<Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false>
292     operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs)
293     {
294       return TriangularProduct
295               <Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false>
296               (lhs.derived(),rhs.m_matrix);
297     }
298 
299     #ifdef EIGEN2_SUPPORT
300     template<typename OtherDerived>
301     struct eigen2_product_return_type
302     {
303       typedef typename TriangularView<MatrixType,Mode>::DenseMatrixType DenseMatrixType;
304       typedef typename OtherDerived::PlainObject::DenseType OtherPlainObject;
305       typedef typename ProductReturnType<DenseMatrixType, OtherPlainObject>::Type ProdRetType;
306       typedef typename ProdRetType::PlainObject type;
307     };
308     template<typename OtherDerived>
309     const typename eigen2_product_return_type<OtherDerived>::type
310     operator*(const EigenBase<OtherDerived>& rhs) const
311     {
312       typename OtherDerived::PlainObject::DenseType rhsPlainObject;
313       rhs.evalTo(rhsPlainObject);
314       return this->toDenseMatrix() * rhsPlainObject;
315     }
316     template<typename OtherMatrixType>
317     bool isApprox(const TriangularView<OtherMatrixType, Mode>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const
318     {
319       return this->toDenseMatrix().isApprox(other.toDenseMatrix(), precision);
320     }
321     template<typename OtherDerived>
322     bool isApprox(const MatrixBase<OtherDerived>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const
323     {
324       return this->toDenseMatrix().isApprox(other, precision);
325     }
326     #endif // EIGEN2_SUPPORT
327 
328     template<int Side, typename Other>
329     inline const internal::triangular_solve_retval<Side,TriangularView, Other>
330     solve(const MatrixBase<Other>& other) const;
331 
332     template<int Side, typename OtherDerived>
333     void solveInPlace(const MatrixBase<OtherDerived>& other) const;
334 
335     template<typename Other>
336     inline const internal::triangular_solve_retval<OnTheLeft,TriangularView, Other>
337     solve(const MatrixBase<Other>& other) const
338     { return solve<OnTheLeft>(other); }
339 
340     template<typename OtherDerived>
341     void solveInPlace(const MatrixBase<OtherDerived>& other) const
342     { return solveInPlace<OnTheLeft>(other); }
343 
344     const SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() const
345     {
346       EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
347       return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
348     }
349     SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView()
350     {
351       EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
352       return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
353     }
354 
355     template<typename OtherDerived>
356     void swap(TriangularBase<OtherDerived> const & other)
357     {
358       TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived());
359     }
360 
361     template<typename OtherDerived>
362     void swap(MatrixBase<OtherDerived> const & other)
363     {
364       SwapWrapper<MatrixType> swaper(const_cast<MatrixType&>(m_matrix));
365       TriangularView<SwapWrapper<MatrixType>,Mode>(swaper).lazyAssign(other.derived());
366     }
367 
368     Scalar determinant() const
369     {
370       if (Mode & UnitDiag)
371         return 1;
372       else if (Mode & ZeroDiag)
373         return 0;
374       else
375         return m_matrix.diagonal().prod();
376     }
377 
378     // TODO simplify the following:
379     template<typename ProductDerived, typename Lhs, typename Rhs>
380     EIGEN_STRONG_INLINE TriangularView& operator=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
381     {
382       setZero();
383       return assignProduct(other.derived(),1);
384     }
385 
386     template<typename ProductDerived, typename Lhs, typename Rhs>
387     EIGEN_STRONG_INLINE TriangularView& operator+=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
388     {
389       return assignProduct(other.derived(),1);
390     }
391 
392     template<typename ProductDerived, typename Lhs, typename Rhs>
393     EIGEN_STRONG_INLINE TriangularView& operator-=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
394     {
395       return assignProduct(other.derived(),-1);
396     }
397 
398 
399     template<typename ProductDerived>
400     EIGEN_STRONG_INLINE TriangularView& operator=(const ScaledProduct<ProductDerived>& other)
401     {
402       setZero();
403       return assignProduct(other.derived(),other.alpha());
404     }
405 
406     template<typename ProductDerived>
407     EIGEN_STRONG_INLINE TriangularView& operator+=(const ScaledProduct<ProductDerived>& other)
408     {
409       return assignProduct(other.derived(),other.alpha());
410     }
411 
412     template<typename ProductDerived>
413     EIGEN_STRONG_INLINE TriangularView& operator-=(const ScaledProduct<ProductDerived>& other)
414     {
415       return assignProduct(other.derived(),-other.alpha());
416     }
417 
418   protected:
419 
420     template<typename ProductDerived, typename Lhs, typename Rhs>
421     EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase<ProductDerived, Lhs,Rhs>& prod, const Scalar& alpha);
422 
423     template<int Mode, bool LhsIsTriangular,
424          typename Lhs, bool LhsIsVector,
425          typename Rhs, bool RhsIsVector>
426     EIGEN_STRONG_INLINE TriangularView& assignProduct(const TriangularProduct<Mode, LhsIsTriangular, Lhs, LhsIsVector, Rhs, RhsIsVector>& prod, const Scalar& alpha)
427     {
428       lazyAssign(alpha*prod.eval());
429       return *this;
430     }
431 
432     MatrixTypeNested m_matrix;
433 };
434 
435 /***************************************************************************
436 * Implementation of triangular evaluation/assignment
437 ***************************************************************************/
438 
439 namespace internal {
440 
441 template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount, bool ClearOpposite>
442 struct triangular_assignment_selector
443 {
444   enum {
445     col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
446     row = (UnrollCount-1) % Derived1::RowsAtCompileTime
447   };
448 
449   typedef typename Derived1::Scalar Scalar;
450 
451   static inline void run(Derived1 &dst, const Derived2 &src)
452   {
453     triangular_assignment_selector<Derived1, Derived2, Mode, UnrollCount-1, ClearOpposite>::run(dst, src);
454 
455     eigen_assert( Mode == Upper || Mode == Lower
456             || Mode == StrictlyUpper || Mode == StrictlyLower
457             || Mode == UnitUpper || Mode == UnitLower);
458     if((Mode == Upper && row <= col)
459     || (Mode == Lower && row >= col)
460     || (Mode == StrictlyUpper && row < col)
461     || (Mode == StrictlyLower && row > col)
462     || (Mode == UnitUpper && row < col)
463     || (Mode == UnitLower && row > col))
464       dst.copyCoeff(row, col, src);
465     else if(ClearOpposite)
466     {
467       if (Mode&UnitDiag && row==col)
468         dst.coeffRef(row, col) = Scalar(1);
469       else
470         dst.coeffRef(row, col) = Scalar(0);
471     }
472   }
473 };
474 
475 // prevent buggy user code from causing an infinite recursion
476 template<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite>
477 struct triangular_assignment_selector<Derived1, Derived2, Mode, 0, ClearOpposite>
478 {
479   static inline void run(Derived1 &, const Derived2 &) {}
480 };
481 
482 template<typename Derived1, typename Derived2, bool ClearOpposite>
483 struct triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearOpposite>
484 {
485   typedef typename Derived1::Index Index;
486   typedef typename Derived1::Scalar Scalar;
487   static inline void run(Derived1 &dst, const Derived2 &src)
488   {
489     for(Index j = 0; j < dst.cols(); ++j)
490     {
491       Index maxi = (std::min)(j, dst.rows()-1);
492       for(Index i = 0; i <= maxi; ++i)
493         dst.copyCoeff(i, j, src);
494       if (ClearOpposite)
495         for(Index i = maxi+1; i < dst.rows(); ++i)
496           dst.coeffRef(i, j) = Scalar(0);
497     }
498   }
499 };
500 
501 template<typename Derived1, typename Derived2, bool ClearOpposite>
502 struct triangular_assignment_selector<Derived1, Derived2, Lower, Dynamic, ClearOpposite>
503 {
504   typedef typename Derived1::Index Index;
505   static inline void run(Derived1 &dst, const Derived2 &src)
506   {
507     for(Index j = 0; j < dst.cols(); ++j)
508     {
509       for(Index i = j; i < dst.rows(); ++i)
510         dst.copyCoeff(i, j, src);
511       Index maxi = (std::min)(j, dst.rows());
512       if (ClearOpposite)
513         for(Index i = 0; i < maxi; ++i)
514           dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
515     }
516   }
517 };
518 
519 template<typename Derived1, typename Derived2, bool ClearOpposite>
520 struct triangular_assignment_selector<Derived1, Derived2, StrictlyUpper, Dynamic, ClearOpposite>
521 {
522   typedef typename Derived1::Index Index;
523   typedef typename Derived1::Scalar Scalar;
524   static inline void run(Derived1 &dst, const Derived2 &src)
525   {
526     for(Index j = 0; j < dst.cols(); ++j)
527     {
528       Index maxi = (std::min)(j, dst.rows());
529       for(Index i = 0; i < maxi; ++i)
530         dst.copyCoeff(i, j, src);
531       if (ClearOpposite)
532         for(Index i = maxi; i < dst.rows(); ++i)
533           dst.coeffRef(i, j) = Scalar(0);
534     }
535   }
536 };
537 
538 template<typename Derived1, typename Derived2, bool ClearOpposite>
539 struct triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic, ClearOpposite>
540 {
541   typedef typename Derived1::Index Index;
542   static inline void run(Derived1 &dst, const Derived2 &src)
543   {
544     for(Index j = 0; j < dst.cols(); ++j)
545     {
546       for(Index i = j+1; i < dst.rows(); ++i)
547         dst.copyCoeff(i, j, src);
548       Index maxi = (std::min)(j, dst.rows()-1);
549       if (ClearOpposite)
550         for(Index i = 0; i <= maxi; ++i)
551           dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
552     }
553   }
554 };
555 
556 template<typename Derived1, typename Derived2, bool ClearOpposite>
557 struct triangular_assignment_selector<Derived1, Derived2, UnitUpper, Dynamic, ClearOpposite>
558 {
559   typedef typename Derived1::Index Index;
560   static inline void run(Derived1 &dst, const Derived2 &src)
561   {
562     for(Index j = 0; j < dst.cols(); ++j)
563     {
564       Index maxi = (std::min)(j, dst.rows());
565       for(Index i = 0; i < maxi; ++i)
566         dst.copyCoeff(i, j, src);
567       if (ClearOpposite)
568       {
569         for(Index i = maxi+1; i < dst.rows(); ++i)
570           dst.coeffRef(i, j) = 0;
571       }
572     }
573     dst.diagonal().setOnes();
574   }
575 };
576 template<typename Derived1, typename Derived2, bool ClearOpposite>
577 struct triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, ClearOpposite>
578 {
579   typedef typename Derived1::Index Index;
580   static inline void run(Derived1 &dst, const Derived2 &src)
581   {
582     for(Index j = 0; j < dst.cols(); ++j)
583     {
584       Index maxi = (std::min)(j, dst.rows());
585       for(Index i = maxi+1; i < dst.rows(); ++i)
586         dst.copyCoeff(i, j, src);
587       if (ClearOpposite)
588       {
589         for(Index i = 0; i < maxi; ++i)
590           dst.coeffRef(i, j) = 0;
591       }
592     }
593     dst.diagonal().setOnes();
594   }
595 };
596 
597 } // end namespace internal
598 
599 // FIXME should we keep that possibility
600 template<typename MatrixType, unsigned int Mode>
601 template<typename OtherDerived>
602 inline TriangularView<MatrixType, Mode>&
603 TriangularView<MatrixType, Mode>::operator=(const MatrixBase<OtherDerived>& other)
604 {
605   if(OtherDerived::Flags & EvalBeforeAssigningBit)
606   {
607     typename internal::plain_matrix_type<OtherDerived>::type other_evaluated(other.rows(), other.cols());
608     other_evaluated.template triangularView<Mode>().lazyAssign(other.derived());
609     lazyAssign(other_evaluated);
610   }
611   else
612     lazyAssign(other.derived());
613   return *this;
614 }
615 
616 // FIXME should we keep that possibility
617 template<typename MatrixType, unsigned int Mode>
618 template<typename OtherDerived>
619 void TriangularView<MatrixType, Mode>::lazyAssign(const MatrixBase<OtherDerived>& other)
620 {
621   enum {
622     unroll = MatrixType::SizeAtCompileTime != Dynamic
623           && internal::traits<OtherDerived>::CoeffReadCost != Dynamic
624           && MatrixType::SizeAtCompileTime*internal::traits<OtherDerived>::CoeffReadCost/2 <= EIGEN_UNROLLING_LIMIT
625   };
626   eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
627 
628   internal::triangular_assignment_selector
629     <MatrixType, OtherDerived, int(Mode),
630     unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic,
631     false // do not change the opposite triangular part
632     >::run(m_matrix.const_cast_derived(), other.derived());
633 }
634 
635 
636 
637 template<typename MatrixType, unsigned int Mode>
638 template<typename OtherDerived>
639 inline TriangularView<MatrixType, Mode>&
640 TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>& other)
641 {
642   eigen_assert(Mode == int(OtherDerived::Mode));
643   if(internal::traits<OtherDerived>::Flags & EvalBeforeAssigningBit)
644   {
645     typename OtherDerived::DenseMatrixType other_evaluated(other.rows(), other.cols());
646     other_evaluated.template triangularView<Mode>().lazyAssign(other.derived().nestedExpression());
647     lazyAssign(other_evaluated);
648   }
649   else
650     lazyAssign(other.derived().nestedExpression());
651   return *this;
652 }
653 
654 template<typename MatrixType, unsigned int Mode>
655 template<typename OtherDerived>
656 void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDerived>& other)
657 {
658   enum {
659     unroll = MatrixType::SizeAtCompileTime != Dynamic
660                    && internal::traits<OtherDerived>::CoeffReadCost != Dynamic
661                    && MatrixType::SizeAtCompileTime * internal::traits<OtherDerived>::CoeffReadCost / 2
662                         <= EIGEN_UNROLLING_LIMIT
663   };
664   eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
665 
666   internal::triangular_assignment_selector
667     <MatrixType, OtherDerived, int(Mode),
668     unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic,
669     false // preserve the opposite triangular part
670     >::run(m_matrix.const_cast_derived(), other.derived().nestedExpression());
671 }
672 
673 /***************************************************************************
674 * Implementation of TriangularBase methods
675 ***************************************************************************/
676 
677 /** Assigns a triangular or selfadjoint matrix to a dense matrix.
678   * If the matrix is triangular, the opposite part is set to zero. */
679 template<typename Derived>
680 template<typename DenseDerived>
681 void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
682 {
683   if(internal::traits<Derived>::Flags & EvalBeforeAssigningBit)
684   {
685     typename internal::plain_matrix_type<Derived>::type other_evaluated(rows(), cols());
686     evalToLazy(other_evaluated);
687     other.derived().swap(other_evaluated);
688   }
689   else
690     evalToLazy(other.derived());
691 }
692 
693 /** Assigns a triangular or selfadjoint matrix to a dense matrix.
694   * If the matrix is triangular, the opposite part is set to zero. */
695 template<typename Derived>
696 template<typename DenseDerived>
697 void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
698 {
699   enum {
700     unroll = DenseDerived::SizeAtCompileTime != Dynamic
701                    && internal::traits<Derived>::CoeffReadCost != Dynamic
702                    && DenseDerived::SizeAtCompileTime * internal::traits<Derived>::CoeffReadCost / 2
703                         <= EIGEN_UNROLLING_LIMIT
704   };
705   other.derived().resize(this->rows(), this->cols());
706 
707   internal::triangular_assignment_selector
708     <DenseDerived, typename internal::traits<Derived>::MatrixTypeNestedCleaned, Derived::Mode,
709     unroll ? int(DenseDerived::SizeAtCompileTime) : Dynamic,
710     true // clear the opposite triangular part
711     >::run(other.derived(), derived().nestedExpression());
712 }
713 
714 /***************************************************************************
715 * Implementation of TriangularView methods
716 ***************************************************************************/
717 
718 /***************************************************************************
719 * Implementation of MatrixBase methods
720 ***************************************************************************/
721 
722 #ifdef EIGEN2_SUPPORT
723 
724 // implementation of part<>(), including the SelfAdjoint case.
725 
726 namespace internal {
727 template<typename MatrixType, unsigned int Mode>
728 struct eigen2_part_return_type
729 {
730   typedef TriangularView<MatrixType, Mode> type;
731 };
732 
733 template<typename MatrixType>
734 struct eigen2_part_return_type<MatrixType, SelfAdjoint>
735 {
736   typedef SelfAdjointView<MatrixType, Upper> type;
737 };
738 }
739 
740 /** \deprecated use MatrixBase::triangularView() */
741 template<typename Derived>
742 template<unsigned int Mode>
743 const typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part() const
744 {
745   return derived();
746 }
747 
748 /** \deprecated use MatrixBase::triangularView() */
749 template<typename Derived>
750 template<unsigned int Mode>
751 typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part()
752 {
753   return derived();
754 }
755 #endif
756 
757 /**
758   * \returns an expression of a triangular view extracted from the current matrix
759   *
760   * The parameter \a Mode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper,
761   * \c #Lower, \c #StrictlyLower, \c #UnitLower.
762   *
763   * Example: \include MatrixBase_extract.cpp
764   * Output: \verbinclude MatrixBase_extract.out
765   *
766   * \sa class TriangularView
767   */
768 template<typename Derived>
769 template<unsigned int Mode>
770 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
771 MatrixBase<Derived>::triangularView()
772 {
773   return derived();
774 }
775 
776 /** This is the const version of MatrixBase::triangularView() */
777 template<typename Derived>
778 template<unsigned int Mode>
779 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
780 MatrixBase<Derived>::triangularView() const
781 {
782   return derived();
783 }
784 
785 /** \returns true if *this is approximately equal to an upper triangular matrix,
786   *          within the precision given by \a prec.
787   *
788   * \sa isLowerTriangular()
789   */
790 template<typename Derived>
791 bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const
792 {
793   using std::abs;
794   RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
795   for(Index j = 0; j < cols(); ++j)
796   {
797     Index maxi = (std::min)(j, rows()-1);
798     for(Index i = 0; i <= maxi; ++i)
799     {
800       RealScalar absValue = abs(coeff(i,j));
801       if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
802     }
803   }
804   RealScalar threshold = maxAbsOnUpperPart * prec;
805   for(Index j = 0; j < cols(); ++j)
806     for(Index i = j+1; i < rows(); ++i)
807       if(abs(coeff(i, j)) > threshold) return false;
808   return true;
809 }
810 
811 /** \returns true if *this is approximately equal to a lower triangular matrix,
812   *          within the precision given by \a prec.
813   *
814   * \sa isUpperTriangular()
815   */
816 template<typename Derived>
817 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const
818 {
819   using std::abs;
820   RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
821   for(Index j = 0; j < cols(); ++j)
822     for(Index i = j; i < rows(); ++i)
823     {
824       RealScalar absValue = abs(coeff(i,j));
825       if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
826     }
827   RealScalar threshold = maxAbsOnLowerPart * prec;
828   for(Index j = 1; j < cols(); ++j)
829   {
830     Index maxi = (std::min)(j, rows()-1);
831     for(Index i = 0; i < maxi; ++i)
832       if(abs(coeff(i, j)) > threshold) return false;
833   }
834   return true;
835 }
836 
837 } // end namespace Eigen
838 
839 #endif // EIGEN_TRIANGULARMATRIX_H
840