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