Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
TriangularMatrix.h
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
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),
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
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
113 template<typename DenseDerived>
114 EIGEN_DEVICE_FUNC
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
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;
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
225 EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
226 inline Index rows() const EIGEN_NOEXCEPT { return m_matrix.rows(); }
228 EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
229 inline Index cols() const EIGEN_NOEXCEPT { return m_matrix.cols(); }
230
232 EIGEN_DEVICE_FUNC
233 const NestedExpression& nestedExpression() const { return m_matrix; }
234
236 EIGEN_DEVICE_FUNC
237 NestedExpression& nestedExpression() { return m_matrix; }
238
241 EIGEN_DEVICE_FUNC
242 inline const ConjugateReturnType conjugate() const
243 { return ConjugateReturnType(m_matrix.conjugate()); }
244
248 template<bool Cond>
249 EIGEN_DEVICE_FUNC
250 inline typename internal::conditional<Cond,ConjugateReturnType,ConstTriangularView>::type
252 {
253 typedef typename internal::conditional<Cond,ConjugateReturnType,ConstTriangularView>::type ReturnType;
254 return ReturnType(m_matrix.template conjugateIf<Cond>());
255 }
256
259 EIGEN_DEVICE_FUNC
260 inline const AdjointReturnType adjoint() const
261 { return AdjointReturnType(m_matrix.adjoint()); }
262
265 EIGEN_DEVICE_FUNC
267 {
268 EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
269 typename MatrixType::TransposeReturnType tmp(m_matrix);
270 return TransposeReturnType(tmp);
271 }
272
275 EIGEN_DEVICE_FUNC
278 return ConstTransposeReturnType(m_matrix.transpose());
279 }
280
281 template<typename Other>
282 EIGEN_DEVICE_FUNC
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
302 EIGEN_DEVICE_FUNC
304 {
305 EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
307 }
308
310 EIGEN_DEVICE_FUNC
312 {
313 EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
315 }
316
317
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
345template<typename _MatrixType, unsigned int _Mode> class TriangularViewImpl<_MatrixType,_Mode,Dense>
346 : public TriangularBase<TriangularView<_MatrixType, _Mode> >
347{
348 public:
349
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
371 EIGEN_DEVICE_FUNC
372 inline Index outerStride() const { return derived().nestedExpression().outerStride(); }
375 EIGEN_DEVICE_FUNC
376 inline Index innerStride() const { return derived().nestedExpression().innerStride(); }
377
379 template<typename Other>
380 EIGEN_DEVICE_FUNC
382 internal::call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op<Scalar,typename Other::Scalar>());
383 return derived();
384 }
386 template<typename Other>
387 EIGEN_DEVICE_FUNC
389 internal::call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op<Scalar,typename Other::Scalar>());
390 return derived();
391 }
392
394 EIGEN_DEVICE_FUNC
395 TriangularViewType& operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() * other; }
397 EIGEN_DEVICE_FUNC
398 TriangularViewType& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() / other; }
399
401 EIGEN_DEVICE_FUNC
402 void fill(const Scalar& value) { setConstant(value); }
404 EIGEN_DEVICE_FUNC
405 TriangularViewType& setConstant(const Scalar& value)
406 { return *this = MatrixType::Constant(derived().rows(), derived().cols(), value); }
408 EIGEN_DEVICE_FUNC
409 TriangularViewType& setZero() { return setConstant(Scalar(0)); }
411 EIGEN_DEVICE_FUNC
412 TriangularViewType& setOnes() { return setConstant(Scalar(1)); }
413
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
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
436 template<typename OtherDerived>
437 EIGEN_DEVICE_FUNC
439
441 template<typename OtherDerived>
442 EIGEN_DEVICE_FUNC
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>
452 EIGEN_DEPRECATED EIGEN_DEVICE_FUNC
453 void lazyAssign(const TriangularBase<OtherDerived>& other);
454
455 template<typename OtherDerived>
457 EIGEN_DEPRECATED EIGEN_DEVICE_FUNC
458 void lazyAssign(const MatrixBase<OtherDerived>& other);
459#endif
460
462 template<typename OtherDerived>
463 EIGEN_DEVICE_FUNC
466 {
467 return Product<TriangularViewType,OtherDerived>(derived(), rhs.derived());
468 }
469
471 template<typename OtherDerived> friend
472 EIGEN_DEVICE_FUNC
474 operator*(const MatrixBase<OtherDerived>& lhs, const TriangularViewImpl& rhs)
475 {
476 return Product<OtherDerived,TriangularViewType>(lhs.derived(),rhs.derived());
477 }
478
502 template<int Side, typename Other>
503 inline const internal::triangular_solve_retval<Side,TriangularViewType, Other>
504 solve(const MatrixBase<Other>& other) const;
505
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
525 template<typename OtherDerived>
526 EIGEN_DEVICE_FUNC
527#ifdef EIGEN_PARSED_BY_DOXYGEN
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
538 template<typename OtherDerived>
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
614template<typename Derived>
615template<typename DenseDerived>
617{
618 evalToLazy(other.derived());
619}
620
621/***************************************************************************
622* Implementation of TriangularView methods
623***************************************************************************/
624
625/***************************************************************************
626* Implementation of MatrixBase methods
627***************************************************************************/
628
640template<typename Derived>
641template<unsigned int Mode>
642EIGEN_DEVICE_FUNC
643typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
645{
646 return typename TriangularViewReturnType<Mode>::Type(derived());
647}
648
650template<typename Derived>
651template<unsigned int Mode>
652EIGEN_DEVICE_FUNC
653typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
655{
656 return typename ConstTriangularViewReturnType<Mode>::Type(derived());
657}
658
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
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
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
949template<typename Derived>
950template<typename DenseDerived>
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
Base class for all dense matrices, vectors, and arrays.
Definition: DenseBase.h:47
Derived & derived()
Definition: EigenBase.h:46
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:65
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:75
Expression of a selfadjoint matrix from a triangular part of a dense matrix.
Definition: SelfAdjointView.h:51
Pseudo expression representing a solving operation.
Definition: Solve.h:63
Base class for triangular part in a matrix.
Definition: TriangularMatrix.h:28
void copyCoeff(Index row, Index col, Other &other)
Definition: TriangularMatrix.h:85
void evalTo(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:616
@ SizeAtCompileTime
Definition: TriangularMatrix.h:38
void evalToLazy(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:951
TriangularViewType & setZero()
Definition: TriangularMatrix.h:409
EIGEN_DEPRECATED void swap(MatrixBase< OtherDerived > const &other)
Definition: TriangularMatrix.h:541
TriangularViewType & operator*=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:395
void solveInPlace(const MatrixBase< OtherDerived > &other) const
void swap(TriangularBase< OtherDerived > &other)
Definition: TriangularMatrix.h:528
TriangularViewType & operator=(const TriangularBase< OtherDerived > &other)
Scalar coeff(Index row, Index col) const
Definition: TriangularMatrix.h:418
Scalar & coeffRef(Index row, Index col)
Definition: TriangularMatrix.h:428
const internal::triangular_solve_retval< Side, TriangularViewType, Other > solve(const MatrixBase< Other > &other) const
friend const Product< OtherDerived, TriangularViewType > operator*(const MatrixBase< OtherDerived > &lhs, const TriangularViewImpl &rhs)
Definition: TriangularMatrix.h:474
TriangularViewType & operator=(const MatrixBase< OtherDerived > &other)
TriangularViewType & operator/=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:398
const Product< TriangularViewType, OtherDerived > operator*(const MatrixBase< OtherDerived > &rhs) const
Definition: TriangularMatrix.h:465
void fill(const Scalar &value)
Definition: TriangularMatrix.h:402
TriangularViewType & operator-=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:388
TriangularViewType & setOnes()
Definition: TriangularMatrix.h:412
Index innerStride() const
Definition: TriangularMatrix.h:376
TriangularViewType & operator+=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:381
TriangularViewType & setConstant(const Scalar &value)
Definition: TriangularMatrix.h:405
Index outerStride() const
Definition: TriangularMatrix.h:372
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:189
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: TriangularMatrix.h:229
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: TriangularMatrix.h:226
const SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView() const
Definition: TriangularMatrix.h:311
const ConjugateReturnType conjugate() const
Definition: TriangularMatrix.h:242
NestedExpression & nestedExpression()
Definition: TriangularMatrix.h:237
internal::conditional< Cond, ConjugateReturnType, ConstTriangularView >::type conjugateIf() const
Definition: TriangularMatrix.h:251
const ConstTransposeReturnType transpose() const
Definition: TriangularMatrix.h:276
const AdjointReturnType adjoint() const
Definition: TriangularMatrix.h:260
Scalar determinant() const
Definition: TriangularMatrix.h:321
SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView()
Definition: TriangularMatrix.h:303
TransposeReturnType transpose()
Definition: TriangularMatrix.h:266
const NestedExpression & nestedExpression() const
Definition: TriangularMatrix.h:233
@ StrictlyLower
Definition: Constants.h:221
@ UnitDiag
Definition: Constants.h:213
@ StrictlyUpper
Definition: Constants.h:223
@ UnitLower
Definition: Constants.h:217
@ ZeroDiag
Definition: Constants.h:215
@ SelfAdjoint
Definition: Constants.h:225
@ UnitUpper
Definition: Constants.h:219
@ Lower
Definition: Constants.h:209
@ Upper
Definition: Constants.h:211
const unsigned int LvalueBit
Definition: Constants.h:144
Namespace containing all symbols from the Eigen library.
Definition: Core:141
const int HugeCost
Definition: Constants.h:44
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const int Dynamic
Definition: Constants.h:22
Definition: Constants.h:507
Definition: EigenBase.h:30
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
Derived & derived()
Definition: EigenBase.h:46