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
14namespace Eigen {
15
16namespace internal {
17
18template<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  */
29template<typename Derived> class TriangularBase : public EigenBase<Derived>
30{
31  public:
32
33    enum {
34      Mode = internal::traits<Derived>::Mode,
35      CoeffReadCost = internal::traits<Derived>::CoeffReadCost,
36      RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
37      ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
38      MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
39      MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime
40    };
41    typedef typename internal::traits<Derived>::Scalar Scalar;
42    typedef typename internal::traits<Derived>::StorageKind StorageKind;
43    typedef typename internal::traits<Derived>::Index Index;
44    typedef typename internal::traits<Derived>::DenseMatrixType DenseMatrixType;
45    typedef DenseMatrixType DenseType;
46
47    inline TriangularBase() { eigen_assert(!((Mode&UnitDiag) && (Mode&ZeroDiag))); }
48
49    inline Index rows() const { return derived().rows(); }
50    inline Index cols() const { return derived().cols(); }
51    inline Index outerStride() const { return derived().outerStride(); }
52    inline Index innerStride() const { return derived().innerStride(); }
53
54    inline Scalar coeff(Index row, Index col) const  { return derived().coeff(row,col); }
55    inline Scalar& coeffRef(Index row, Index col) { return derived().coeffRef(row,col); }
56
57    /** \see MatrixBase::copyCoeff(row,col)
58      */
59    template<typename Other>
60    EIGEN_STRONG_INLINE void copyCoeff(Index row, Index col, Other& other)
61    {
62      derived().coeffRef(row, col) = other.coeff(row, col);
63    }
64
65    inline Scalar operator()(Index row, Index col) const
66    {
67      check_coordinates(row, col);
68      return coeff(row,col);
69    }
70    inline Scalar& operator()(Index row, Index col)
71    {
72      check_coordinates(row, col);
73      return coeffRef(row,col);
74    }
75
76    #ifndef EIGEN_PARSED_BY_DOXYGEN
77    inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
78    inline Derived& derived() { return *static_cast<Derived*>(this); }
79    #endif // not EIGEN_PARSED_BY_DOXYGEN
80
81    template<typename DenseDerived>
82    void evalTo(MatrixBase<DenseDerived> &other) const;
83    template<typename DenseDerived>
84    void evalToLazy(MatrixBase<DenseDerived> &other) const;
85
86    DenseMatrixType toDenseMatrix() const
87    {
88      DenseMatrixType res(rows(), cols());
89      evalToLazy(res);
90      return res;
91    }
92
93  protected:
94
95    void check_coordinates(Index row, Index col) const
96    {
97      EIGEN_ONLY_USED_FOR_DEBUG(row);
98      EIGEN_ONLY_USED_FOR_DEBUG(col);
99      eigen_assert(col>=0 && col<cols() && row>=0 && row<rows());
100      const int mode = int(Mode) & ~SelfAdjoint;
101      EIGEN_ONLY_USED_FOR_DEBUG(mode);
102      eigen_assert((mode==Upper && col>=row)
103                || (mode==Lower && col<=row)
104                || ((mode==StrictlyUpper || mode==UnitUpper) && col>row)
105                || ((mode==StrictlyLower || mode==UnitLower) && col<row));
106    }
107
108    #ifdef EIGEN_INTERNAL_DEBUGGING
109    void check_coordinates_internal(Index row, Index col) const
110    {
111      check_coordinates(row, col);
112    }
113    #else
114    void check_coordinates_internal(Index , Index ) const {}
115    #endif
116
117};
118
119/** \class TriangularView
120  * \ingroup Core_Module
121  *
122  * \brief Base class for triangular part in a matrix
123  *
124  * \param MatrixType the type of the object in which we are taking the triangular part
125  * \param Mode the kind of triangular matrix expression to construct. Can be #Upper,
126  *             #Lower, #UnitUpper, #UnitLower, #StrictlyUpper, or #StrictlyLower.
127  *             This is in fact a bit field; it must have either #Upper or #Lower,
128  *             and additionnaly it may have #UnitDiag or #ZeroDiag or neither.
129  *
130  * This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular
131  * matrices one should speak of "trapezoid" parts. This class is the return type
132  * of MatrixBase::triangularView() and most of the time this is the only way it is used.
133  *
134  * \sa MatrixBase::triangularView()
135  */
136namespace internal {
137template<typename MatrixType, unsigned int _Mode>
138struct 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
153template<int Mode, bool LhsIsTriangular,
154         typename Lhs, bool LhsIsVector,
155         typename Rhs, bool RhsIsVector>
156struct TriangularProduct;
157
158template<typename _MatrixType, unsigned int _Mode> class TriangularView
159  : public TriangularBase<TriangularView<_MatrixType, _Mode> >
160{
161  public:
162
163    typedef TriangularBase<TriangularView> Base;
164    typedef typename internal::traits<TriangularView>::Scalar Scalar;
165
166    typedef _MatrixType MatrixType;
167    typedef typename internal::traits<TriangularView>::DenseMatrixType DenseMatrixType;
168    typedef DenseMatrixType PlainObject;
169
170  protected:
171    typedef typename internal::traits<TriangularView>::MatrixTypeNested MatrixTypeNested;
172    typedef typename internal::traits<TriangularView>::MatrixTypeNestedNonRef MatrixTypeNestedNonRef;
173    typedef typename internal::traits<TriangularView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned;
174
175    typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType;
176
177  public:
178    using Base::evalToLazy;
179
180
181    typedef typename internal::traits<TriangularView>::StorageKind StorageKind;
182    typedef typename internal::traits<TriangularView>::Index Index;
183
184    enum {
185      Mode = _Mode,
186      TransposeMode = (Mode & Upper ? Lower : 0)
187                    | (Mode & Lower ? Upper : 0)
188                    | (Mode & (UnitDiag))
189                    | (Mode & (ZeroDiag))
190    };
191
192    inline TriangularView(const MatrixType& matrix) : m_matrix(matrix)
193    {}
194
195    inline Index rows() const { return m_matrix.rows(); }
196    inline Index cols() const { return m_matrix.cols(); }
197    inline Index outerStride() const { return m_matrix.outerStride(); }
198    inline Index innerStride() const { return m_matrix.innerStride(); }
199
200    /** \sa MatrixBase::operator+=() */
201    template<typename Other> TriangularView&  operator+=(const DenseBase<Other>& other) { return *this = m_matrix + other.derived(); }
202    /** \sa MatrixBase::operator-=() */
203    template<typename Other> TriangularView&  operator-=(const DenseBase<Other>& other) { return *this = m_matrix - other.derived(); }
204    /** \sa MatrixBase::operator*=() */
205    TriangularView&  operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix * other; }
206    /** \sa MatrixBase::operator/=() */
207    TriangularView&  operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix / other; }
208
209    /** \sa MatrixBase::fill() */
210    void fill(const Scalar& value) { setConstant(value); }
211    /** \sa MatrixBase::setConstant() */
212    TriangularView& setConstant(const Scalar& value)
213    { return *this = MatrixType::Constant(rows(), cols(), value); }
214    /** \sa MatrixBase::setZero() */
215    TriangularView& setZero() { return setConstant(Scalar(0)); }
216    /** \sa MatrixBase::setOnes() */
217    TriangularView& setOnes() { return setConstant(Scalar(1)); }
218
219    /** \sa MatrixBase::coeff()
220      * \warning the coordinates must fit into the referenced triangular part
221      */
222    inline Scalar coeff(Index row, Index col) const
223    {
224      Base::check_coordinates_internal(row, col);
225      return m_matrix.coeff(row, col);
226    }
227
228    /** \sa MatrixBase::coeffRef()
229      * \warning the coordinates must fit into the referenced triangular part
230      */
231    inline Scalar& coeffRef(Index row, Index col)
232    {
233      Base::check_coordinates_internal(row, col);
234      return m_matrix.const_cast_derived().coeffRef(row, col);
235    }
236
237    const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; }
238    MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); }
239
240    /** Assigns a triangular matrix to a triangular part of a dense matrix */
241    template<typename OtherDerived>
242    TriangularView& operator=(const TriangularBase<OtherDerived>& other);
243
244    template<typename OtherDerived>
245    TriangularView& operator=(const MatrixBase<OtherDerived>& other);
246
247    TriangularView& operator=(const TriangularView& other)
248    { return *this = other.nestedExpression(); }
249
250    template<typename OtherDerived>
251    void lazyAssign(const TriangularBase<OtherDerived>& other);
252
253    template<typename OtherDerived>
254    void lazyAssign(const MatrixBase<OtherDerived>& other);
255
256    /** \sa MatrixBase::conjugate() */
257    inline TriangularView<MatrixConjugateReturnType,Mode> conjugate()
258    { return m_matrix.conjugate(); }
259    /** \sa MatrixBase::conjugate() const */
260    inline const TriangularView<MatrixConjugateReturnType,Mode> conjugate() const
261    { return m_matrix.conjugate(); }
262
263    /** \sa MatrixBase::adjoint() const */
264    inline const TriangularView<const typename MatrixType::AdjointReturnType,TransposeMode> adjoint() const
265    { return m_matrix.adjoint(); }
266
267    /** \sa MatrixBase::transpose() */
268    inline TriangularView<Transpose<MatrixType>,TransposeMode> transpose()
269    {
270      EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
271      return m_matrix.const_cast_derived().transpose();
272    }
273    /** \sa MatrixBase::transpose() const */
274    inline const TriangularView<Transpose<MatrixType>,TransposeMode> transpose() const
275    {
276      return m_matrix.transpose();
277    }
278
279    /** Efficient triangular matrix times vector/matrix product */
280    template<typename OtherDerived>
281    TriangularProduct<Mode, true, MatrixType, false, OtherDerived, OtherDerived::ColsAtCompileTime==1>
282    operator*(const MatrixBase<OtherDerived>& rhs) const
283    {
284      return TriangularProduct
285              <Mode, true, MatrixType, false, OtherDerived, OtherDerived::ColsAtCompileTime==1>
286              (m_matrix, rhs.derived());
287    }
288
289    /** Efficient vector/matrix times triangular matrix product */
290    template<typename OtherDerived> friend
291    TriangularProduct<Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false>
292    operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs)
293    {
294      return TriangularProduct
295              <Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false>
296              (lhs.derived(),rhs.m_matrix);
297    }
298
299    #ifdef EIGEN2_SUPPORT
300    template<typename OtherDerived>
301    struct eigen2_product_return_type
302    {
303      typedef typename TriangularView<MatrixType,Mode>::DenseMatrixType DenseMatrixType;
304      typedef typename OtherDerived::PlainObject::DenseType OtherPlainObject;
305      typedef typename ProductReturnType<DenseMatrixType, OtherPlainObject>::Type ProdRetType;
306      typedef typename ProdRetType::PlainObject type;
307    };
308    template<typename OtherDerived>
309    const typename eigen2_product_return_type<OtherDerived>::type
310    operator*(const EigenBase<OtherDerived>& rhs) const
311    {
312      typename OtherDerived::PlainObject::DenseType rhsPlainObject;
313      rhs.evalTo(rhsPlainObject);
314      return this->toDenseMatrix() * rhsPlainObject;
315    }
316    template<typename OtherMatrixType>
317    bool isApprox(const TriangularView<OtherMatrixType, Mode>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const
318    {
319      return this->toDenseMatrix().isApprox(other.toDenseMatrix(), precision);
320    }
321    template<typename OtherDerived>
322    bool isApprox(const MatrixBase<OtherDerived>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const
323    {
324      return this->toDenseMatrix().isApprox(other, precision);
325    }
326    #endif // EIGEN2_SUPPORT
327
328    template<int Side, typename Other>
329    inline const internal::triangular_solve_retval<Side,TriangularView, Other>
330    solve(const MatrixBase<Other>& other) const;
331
332    template<int Side, typename OtherDerived>
333    void solveInPlace(const MatrixBase<OtherDerived>& other) const;
334
335    template<typename Other>
336    inline const internal::triangular_solve_retval<OnTheLeft,TriangularView, Other>
337    solve(const MatrixBase<Other>& other) const
338    { return solve<OnTheLeft>(other); }
339
340    template<typename OtherDerived>
341    void solveInPlace(const MatrixBase<OtherDerived>& other) const
342    { return solveInPlace<OnTheLeft>(other); }
343
344    const SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() const
345    {
346      EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
347      return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
348    }
349    SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView()
350    {
351      EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
352      return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
353    }
354
355    template<typename OtherDerived>
356    void swap(TriangularBase<OtherDerived> const & other)
357    {
358      TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived());
359    }
360
361    template<typename OtherDerived>
362    void swap(MatrixBase<OtherDerived> const & other)
363    {
364      SwapWrapper<MatrixType> swaper(const_cast<MatrixType&>(m_matrix));
365      TriangularView<SwapWrapper<MatrixType>,Mode>(swaper).lazyAssign(other.derived());
366    }
367
368    Scalar determinant() const
369    {
370      if (Mode & UnitDiag)
371        return 1;
372      else if (Mode & ZeroDiag)
373        return 0;
374      else
375        return m_matrix.diagonal().prod();
376    }
377
378    // TODO simplify the following:
379    template<typename ProductDerived, typename Lhs, typename Rhs>
380    EIGEN_STRONG_INLINE TriangularView& operator=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
381    {
382      setZero();
383      return assignProduct(other,1);
384    }
385
386    template<typename ProductDerived, typename Lhs, typename Rhs>
387    EIGEN_STRONG_INLINE TriangularView& operator+=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
388    {
389      return assignProduct(other,1);
390    }
391
392    template<typename ProductDerived, typename Lhs, typename Rhs>
393    EIGEN_STRONG_INLINE TriangularView& operator-=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
394    {
395      return assignProduct(other,-1);
396    }
397
398
399    template<typename ProductDerived>
400    EIGEN_STRONG_INLINE TriangularView& operator=(const ScaledProduct<ProductDerived>& other)
401    {
402      setZero();
403      return assignProduct(other,other.alpha());
404    }
405
406    template<typename ProductDerived>
407    EIGEN_STRONG_INLINE TriangularView& operator+=(const ScaledProduct<ProductDerived>& other)
408    {
409      return assignProduct(other,other.alpha());
410    }
411
412    template<typename ProductDerived>
413    EIGEN_STRONG_INLINE TriangularView& operator-=(const ScaledProduct<ProductDerived>& other)
414    {
415      return assignProduct(other,-other.alpha());
416    }
417
418  protected:
419
420    template<typename ProductDerived, typename Lhs, typename Rhs>
421    EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase<ProductDerived, Lhs,Rhs>& prod, const Scalar& alpha);
422
423    MatrixTypeNested m_matrix;
424};
425
426/***************************************************************************
427* Implementation of triangular evaluation/assignment
428***************************************************************************/
429
430namespace internal {
431
432template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount, bool ClearOpposite>
433struct triangular_assignment_selector
434{
435  enum {
436    col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
437    row = (UnrollCount-1) % Derived1::RowsAtCompileTime
438  };
439
440  typedef typename Derived1::Scalar Scalar;
441
442  static inline void run(Derived1 &dst, const Derived2 &src)
443  {
444    triangular_assignment_selector<Derived1, Derived2, Mode, UnrollCount-1, ClearOpposite>::run(dst, src);
445
446    eigen_assert( Mode == Upper || Mode == Lower
447            || Mode == StrictlyUpper || Mode == StrictlyLower
448            || Mode == UnitUpper || Mode == UnitLower);
449    if((Mode == Upper && row <= col)
450    || (Mode == Lower && row >= col)
451    || (Mode == StrictlyUpper && row < col)
452    || (Mode == StrictlyLower && row > col)
453    || (Mode == UnitUpper && row < col)
454    || (Mode == UnitLower && row > col))
455      dst.copyCoeff(row, col, src);
456    else if(ClearOpposite)
457    {
458      if (Mode&UnitDiag && row==col)
459        dst.coeffRef(row, col) = Scalar(1);
460      else
461        dst.coeffRef(row, col) = Scalar(0);
462    }
463  }
464};
465
466// prevent buggy user code from causing an infinite recursion
467template<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite>
468struct triangular_assignment_selector<Derived1, Derived2, Mode, 0, ClearOpposite>
469{
470  static inline void run(Derived1 &, const Derived2 &) {}
471};
472
473template<typename Derived1, typename Derived2, bool ClearOpposite>
474struct triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearOpposite>
475{
476  typedef typename Derived1::Index Index;
477  typedef typename Derived1::Scalar Scalar;
478  static inline void run(Derived1 &dst, const Derived2 &src)
479  {
480    for(Index j = 0; j < dst.cols(); ++j)
481    {
482      Index maxi = (std::min)(j, dst.rows()-1);
483      for(Index i = 0; i <= maxi; ++i)
484        dst.copyCoeff(i, j, src);
485      if (ClearOpposite)
486        for(Index i = maxi+1; i < dst.rows(); ++i)
487          dst.coeffRef(i, j) = Scalar(0);
488    }
489  }
490};
491
492template<typename Derived1, typename Derived2, bool ClearOpposite>
493struct triangular_assignment_selector<Derived1, Derived2, Lower, Dynamic, ClearOpposite>
494{
495  typedef typename Derived1::Index Index;
496  static inline void run(Derived1 &dst, const Derived2 &src)
497  {
498    for(Index j = 0; j < dst.cols(); ++j)
499    {
500      for(Index i = j; i < dst.rows(); ++i)
501        dst.copyCoeff(i, j, src);
502      Index maxi = (std::min)(j, dst.rows());
503      if (ClearOpposite)
504        for(Index i = 0; i < maxi; ++i)
505          dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
506    }
507  }
508};
509
510template<typename Derived1, typename Derived2, bool ClearOpposite>
511struct triangular_assignment_selector<Derived1, Derived2, StrictlyUpper, Dynamic, ClearOpposite>
512{
513  typedef typename Derived1::Index Index;
514  typedef typename Derived1::Scalar Scalar;
515  static inline void run(Derived1 &dst, const Derived2 &src)
516  {
517    for(Index j = 0; j < dst.cols(); ++j)
518    {
519      Index maxi = (std::min)(j, dst.rows());
520      for(Index i = 0; i < maxi; ++i)
521        dst.copyCoeff(i, j, src);
522      if (ClearOpposite)
523        for(Index i = maxi; i < dst.rows(); ++i)
524          dst.coeffRef(i, j) = Scalar(0);
525    }
526  }
527};
528
529template<typename Derived1, typename Derived2, bool ClearOpposite>
530struct triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic, ClearOpposite>
531{
532  typedef typename Derived1::Index Index;
533  static inline void run(Derived1 &dst, const Derived2 &src)
534  {
535    for(Index j = 0; j < dst.cols(); ++j)
536    {
537      for(Index i = j+1; i < dst.rows(); ++i)
538        dst.copyCoeff(i, j, src);
539      Index maxi = (std::min)(j, dst.rows()-1);
540      if (ClearOpposite)
541        for(Index i = 0; i <= maxi; ++i)
542          dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
543    }
544  }
545};
546
547template<typename Derived1, typename Derived2, bool ClearOpposite>
548struct triangular_assignment_selector<Derived1, Derived2, UnitUpper, Dynamic, ClearOpposite>
549{
550  typedef typename Derived1::Index Index;
551  static inline void run(Derived1 &dst, const Derived2 &src)
552  {
553    for(Index j = 0; j < dst.cols(); ++j)
554    {
555      Index maxi = (std::min)(j, dst.rows());
556      for(Index i = 0; i < maxi; ++i)
557        dst.copyCoeff(i, j, src);
558      if (ClearOpposite)
559      {
560        for(Index i = maxi+1; i < dst.rows(); ++i)
561          dst.coeffRef(i, j) = 0;
562      }
563    }
564    dst.diagonal().setOnes();
565  }
566};
567template<typename Derived1, typename Derived2, bool ClearOpposite>
568struct triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, ClearOpposite>
569{
570  typedef typename Derived1::Index Index;
571  static inline void run(Derived1 &dst, const Derived2 &src)
572  {
573    for(Index j = 0; j < dst.cols(); ++j)
574    {
575      Index maxi = (std::min)(j, dst.rows());
576      for(Index i = maxi+1; i < dst.rows(); ++i)
577        dst.copyCoeff(i, j, src);
578      if (ClearOpposite)
579      {
580        for(Index i = 0; i < maxi; ++i)
581          dst.coeffRef(i, j) = 0;
582      }
583    }
584    dst.diagonal().setOnes();
585  }
586};
587
588} // end namespace internal
589
590// FIXME should we keep that possibility
591template<typename MatrixType, unsigned int Mode>
592template<typename OtherDerived>
593inline TriangularView<MatrixType, Mode>&
594TriangularView<MatrixType, Mode>::operator=(const MatrixBase<OtherDerived>& other)
595{
596  if(OtherDerived::Flags & EvalBeforeAssigningBit)
597  {
598    typename internal::plain_matrix_type<OtherDerived>::type other_evaluated(other.rows(), other.cols());
599    other_evaluated.template triangularView<Mode>().lazyAssign(other.derived());
600    lazyAssign(other_evaluated);
601  }
602  else
603    lazyAssign(other.derived());
604  return *this;
605}
606
607// FIXME should we keep that possibility
608template<typename MatrixType, unsigned int Mode>
609template<typename OtherDerived>
610void TriangularView<MatrixType, Mode>::lazyAssign(const MatrixBase<OtherDerived>& other)
611{
612  enum {
613    unroll = MatrixType::SizeAtCompileTime != Dynamic
614          && internal::traits<OtherDerived>::CoeffReadCost != Dynamic
615          && MatrixType::SizeAtCompileTime*internal::traits<OtherDerived>::CoeffReadCost/2 <= EIGEN_UNROLLING_LIMIT
616  };
617  eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
618
619  internal::triangular_assignment_selector
620    <MatrixType, OtherDerived, int(Mode),
621    unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic,
622    false // do not change the opposite triangular part
623    >::run(m_matrix.const_cast_derived(), other.derived());
624}
625
626
627
628template<typename MatrixType, unsigned int Mode>
629template<typename OtherDerived>
630inline TriangularView<MatrixType, Mode>&
631TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>& other)
632{
633  eigen_assert(Mode == int(OtherDerived::Mode));
634  if(internal::traits<OtherDerived>::Flags & EvalBeforeAssigningBit)
635  {
636    typename OtherDerived::DenseMatrixType other_evaluated(other.rows(), other.cols());
637    other_evaluated.template triangularView<Mode>().lazyAssign(other.derived().nestedExpression());
638    lazyAssign(other_evaluated);
639  }
640  else
641    lazyAssign(other.derived().nestedExpression());
642  return *this;
643}
644
645template<typename MatrixType, unsigned int Mode>
646template<typename OtherDerived>
647void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDerived>& other)
648{
649  enum {
650    unroll = MatrixType::SizeAtCompileTime != Dynamic
651                   && internal::traits<OtherDerived>::CoeffReadCost != Dynamic
652                   && MatrixType::SizeAtCompileTime * internal::traits<OtherDerived>::CoeffReadCost / 2
653                        <= EIGEN_UNROLLING_LIMIT
654  };
655  eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
656
657  internal::triangular_assignment_selector
658    <MatrixType, OtherDerived, int(Mode),
659    unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic,
660    false // preserve the opposite triangular part
661    >::run(m_matrix.const_cast_derived(), other.derived().nestedExpression());
662}
663
664/***************************************************************************
665* Implementation of TriangularBase methods
666***************************************************************************/
667
668/** Assigns a triangular or selfadjoint matrix to a dense matrix.
669  * If the matrix is triangular, the opposite part is set to zero. */
670template<typename Derived>
671template<typename DenseDerived>
672void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
673{
674  if(internal::traits<Derived>::Flags & EvalBeforeAssigningBit)
675  {
676    typename internal::plain_matrix_type<Derived>::type other_evaluated(rows(), cols());
677    evalToLazy(other_evaluated);
678    other.derived().swap(other_evaluated);
679  }
680  else
681    evalToLazy(other.derived());
682}
683
684/** Assigns a triangular or selfadjoint matrix to a dense matrix.
685  * If the matrix is triangular, the opposite part is set to zero. */
686template<typename Derived>
687template<typename DenseDerived>
688void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
689{
690  enum {
691    unroll = DenseDerived::SizeAtCompileTime != Dynamic
692                   && internal::traits<Derived>::CoeffReadCost != Dynamic
693                   && DenseDerived::SizeAtCompileTime * internal::traits<Derived>::CoeffReadCost / 2
694                        <= EIGEN_UNROLLING_LIMIT
695  };
696  other.derived().resize(this->rows(), this->cols());
697
698  internal::triangular_assignment_selector
699    <DenseDerived, typename internal::traits<Derived>::MatrixTypeNestedCleaned, Derived::Mode,
700    unroll ? int(DenseDerived::SizeAtCompileTime) : Dynamic,
701    true // clear the opposite triangular part
702    >::run(other.derived(), derived().nestedExpression());
703}
704
705/***************************************************************************
706* Implementation of TriangularView methods
707***************************************************************************/
708
709/***************************************************************************
710* Implementation of MatrixBase methods
711***************************************************************************/
712
713#ifdef EIGEN2_SUPPORT
714
715// implementation of part<>(), including the SelfAdjoint case.
716
717namespace internal {
718template<typename MatrixType, unsigned int Mode>
719struct eigen2_part_return_type
720{
721  typedef TriangularView<MatrixType, Mode> type;
722};
723
724template<typename MatrixType>
725struct eigen2_part_return_type<MatrixType, SelfAdjoint>
726{
727  typedef SelfAdjointView<MatrixType, Upper> type;
728};
729}
730
731/** \deprecated use MatrixBase::triangularView() */
732template<typename Derived>
733template<unsigned int Mode>
734const typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part() const
735{
736  return derived();
737}
738
739/** \deprecated use MatrixBase::triangularView() */
740template<typename Derived>
741template<unsigned int Mode>
742typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part()
743{
744  return derived();
745}
746#endif
747
748/**
749  * \returns an expression of a triangular view extracted from the current matrix
750  *
751  * The parameter \a Mode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper,
752  * \c #Lower, \c #StrictlyLower, \c #UnitLower.
753  *
754  * Example: \include MatrixBase_extract.cpp
755  * Output: \verbinclude MatrixBase_extract.out
756  *
757  * \sa class TriangularView
758  */
759template<typename Derived>
760template<unsigned int Mode>
761typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
762MatrixBase<Derived>::triangularView()
763{
764  return derived();
765}
766
767/** This is the const version of MatrixBase::triangularView() */
768template<typename Derived>
769template<unsigned int Mode>
770typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
771MatrixBase<Derived>::triangularView() const
772{
773  return derived();
774}
775
776/** \returns true if *this is approximately equal to an upper triangular matrix,
777  *          within the precision given by \a prec.
778  *
779  * \sa isLowerTriangular()
780  */
781template<typename Derived>
782bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const
783{
784  using std::abs;
785  RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
786  for(Index j = 0; j < cols(); ++j)
787  {
788    Index maxi = (std::min)(j, rows()-1);
789    for(Index i = 0; i <= maxi; ++i)
790    {
791      RealScalar absValue = abs(coeff(i,j));
792      if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
793    }
794  }
795  RealScalar threshold = maxAbsOnUpperPart * prec;
796  for(Index j = 0; j < cols(); ++j)
797    for(Index i = j+1; i < rows(); ++i)
798      if(abs(coeff(i, j)) > threshold) return false;
799  return true;
800}
801
802/** \returns true if *this is approximately equal to a lower triangular matrix,
803  *          within the precision given by \a prec.
804  *
805  * \sa isUpperTriangular()
806  */
807template<typename Derived>
808bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const
809{
810  using std::abs;
811  RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
812  for(Index j = 0; j < cols(); ++j)
813    for(Index i = j; i < rows(); ++i)
814    {
815      RealScalar absValue = abs(coeff(i,j));
816      if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
817    }
818  RealScalar threshold = maxAbsOnLowerPart * prec;
819  for(Index j = 1; j < cols(); ++j)
820  {
821    Index maxi = (std::min)(j, rows()-1);
822    for(Index i = 0; i < maxi; ++i)
823      if(abs(coeff(i, j)) > threshold) return false;
824  }
825  return true;
826}
827
828} // end namespace Eigen
829
830#endif // EIGEN_TRIANGULARMATRIX_H
831