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