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

source code of qtmultimedia/src/3rdparty/eigen/Eigen/src/Core/TriangularMatrix.h