Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.3.9
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 
14 namespace Eigen {
15 
16 namespace internal {
17 
18 template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval;
19 
20 }
21 
27 template<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(!((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 
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
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 
167 namespace internal {
168 template<typename MatrixType, unsigned int _Mode>
169 struct 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 
184 template<typename _MatrixType, unsigned int _Mode, typename StorageKind> class TriangularViewImpl;
185 
186 template<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  EIGEN_INHERIT_ASSIGNMENT_OPERATORS(TriangularView)
221 
222 
223  EIGEN_DEVICE_FUNC
224  inline Index rows() const { return m_matrix.rows(); }
226  EIGEN_DEVICE_FUNC
227  inline Index cols() const { return m_matrix.cols(); }
228 
230  EIGEN_DEVICE_FUNC
231  const NestedExpression& nestedExpression() const { return m_matrix; }
232 
234  EIGEN_DEVICE_FUNC
235  NestedExpression& nestedExpression() { return m_matrix; }
236 
237  typedef TriangularView<const MatrixConjugateReturnType,Mode> ConjugateReturnType;
239  EIGEN_DEVICE_FUNC
240  inline const ConjugateReturnType conjugate() const
241  { return ConjugateReturnType(m_matrix.conjugate()); }
242 
245  EIGEN_DEVICE_FUNC
246  inline const AdjointReturnType adjoint() const
247  { return AdjointReturnType(m_matrix.adjoint()); }
248 
251  EIGEN_DEVICE_FUNC
253  {
254  EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
255  typename MatrixType::TransposeReturnType tmp(m_matrix);
256  return TransposeReturnType(tmp);
257  }
258 
261  EIGEN_DEVICE_FUNC
262  inline const ConstTransposeReturnType transpose() const
263  {
264  return ConstTransposeReturnType(m_matrix.transpose());
265  }
266 
267  template<typename Other>
268  EIGEN_DEVICE_FUNC
269  inline const Solve<TriangularView, Other>
270  solve(const MatrixBase<Other>& other) const
271  { return Solve<TriangularView, Other>(*this, other.derived()); }
272 
273  // workaround MSVC ICE
274  #if EIGEN_COMP_MSVC
275  template<int Side, typename Other>
276  EIGEN_DEVICE_FUNC
277  inline const internal::triangular_solve_retval<Side,TriangularView, Other>
278  solve(const MatrixBase<Other>& other) const
279  { return Base::template solve<Side>(other); }
280  #else
281  using Base::solve;
282  #endif
283 
288  EIGEN_DEVICE_FUNC
290  {
291  EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
293  }
294 
296  EIGEN_DEVICE_FUNC
298  {
299  EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
301  }
302 
303 
306  EIGEN_DEVICE_FUNC
307  Scalar determinant() const
308  {
309  if (Mode & UnitDiag)
310  return 1;
311  else if (Mode & ZeroDiag)
312  return 0;
313  else
314  return m_matrix.diagonal().prod();
315  }
316 
317  protected:
318 
319  MatrixTypeNested m_matrix;
320 };
321 
331 template<typename _MatrixType, unsigned int _Mode> class TriangularViewImpl<_MatrixType,_Mode,Dense>
332  : public TriangularBase<TriangularView<_MatrixType, _Mode> >
333 {
334  public:
335 
338  typedef typename internal::traits<TriangularViewType>::Scalar Scalar;
339 
340  typedef _MatrixType MatrixType;
341  typedef typename MatrixType::PlainObject DenseMatrixType;
342  typedef DenseMatrixType PlainObject;
343 
344  public:
345  using Base::evalToLazy;
346  using Base::derived;
347 
348  typedef typename internal::traits<TriangularViewType>::StorageKind StorageKind;
349 
350  enum {
351  Mode = _Mode,
352  Flags = internal::traits<TriangularViewType>::Flags
353  };
354 
357  EIGEN_DEVICE_FUNC
358  inline Index outerStride() const { return derived().nestedExpression().outerStride(); }
361  EIGEN_DEVICE_FUNC
362  inline Index innerStride() const { return derived().nestedExpression().innerStride(); }
363 
365  template<typename Other>
366  EIGEN_DEVICE_FUNC
368  internal::call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op<Scalar,typename Other::Scalar>());
369  return derived();
370  }
372  template<typename Other>
373  EIGEN_DEVICE_FUNC
375  internal::call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op<Scalar,typename Other::Scalar>());
376  return derived();
377  }
378 
380  EIGEN_DEVICE_FUNC
381  TriangularViewType& operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() * other; }
383  EIGEN_DEVICE_FUNC
384  TriangularViewType& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() / other; }
385 
387  EIGEN_DEVICE_FUNC
388  void fill(const Scalar& value) { setConstant(value); }
390  EIGEN_DEVICE_FUNC
391  TriangularViewType& setConstant(const Scalar& value)
392  { return *this = MatrixType::Constant(derived().rows(), derived().cols(), value); }
394  EIGEN_DEVICE_FUNC
395  TriangularViewType& setZero() { return setConstant(Scalar(0)); }
397  EIGEN_DEVICE_FUNC
398  TriangularViewType& setOnes() { return setConstant(Scalar(1)); }
399 
403  EIGEN_DEVICE_FUNC
404  inline Scalar coeff(Index row, Index col) const
405  {
406  Base::check_coordinates_internal(row, col);
407  return derived().nestedExpression().coeff(row, col);
408  }
409 
413  EIGEN_DEVICE_FUNC
414  inline Scalar& coeffRef(Index row, Index col)
415  {
416  EIGEN_STATIC_ASSERT_LVALUE(TriangularViewType);
417  Base::check_coordinates_internal(row, col);
418  return derived().nestedExpression().coeffRef(row, col);
419  }
420 
422  template<typename OtherDerived>
423  EIGEN_DEVICE_FUNC
425 
427  template<typename OtherDerived>
428  EIGEN_DEVICE_FUNC
430 
431 #ifndef EIGEN_PARSED_BY_DOXYGEN
432  EIGEN_DEVICE_FUNC
433  TriangularViewType& operator=(const TriangularViewImpl& other)
434  { return *this = other.derived().nestedExpression(); }
435 
437  template<typename OtherDerived>
438  EIGEN_DEVICE_FUNC
439  void lazyAssign(const TriangularBase<OtherDerived>& other);
440 
442  template<typename OtherDerived>
443  EIGEN_DEVICE_FUNC
444  void lazyAssign(const MatrixBase<OtherDerived>& other);
445 #endif
446 
448  template<typename OtherDerived>
449  EIGEN_DEVICE_FUNC
452  {
453  return Product<TriangularViewType,OtherDerived>(derived(), rhs.derived());
454  }
455 
457  template<typename OtherDerived> friend
458  EIGEN_DEVICE_FUNC
460  operator*(const MatrixBase<OtherDerived>& lhs, const TriangularViewImpl& rhs)
461  {
462  return Product<OtherDerived,TriangularViewType>(lhs.derived(),rhs.derived());
463  }
464 
488  template<int Side, typename Other>
489  EIGEN_DEVICE_FUNC
490  inline const internal::triangular_solve_retval<Side,TriangularViewType, Other>
491  solve(const MatrixBase<Other>& other) const;
492 
502  template<int Side, typename OtherDerived>
503  EIGEN_DEVICE_FUNC
504  void solveInPlace(const MatrixBase<OtherDerived>& other) const;
505 
506  template<typename OtherDerived>
507  EIGEN_DEVICE_FUNC
508  void solveInPlace(const MatrixBase<OtherDerived>& other) const
509  { return solveInPlace<OnTheLeft>(other); }
510 
512  template<typename OtherDerived>
513  EIGEN_DEVICE_FUNC
514 #ifdef EIGEN_PARSED_BY_DOXYGEN
516 #else
517  void swap(TriangularBase<OtherDerived> const & other)
518 #endif
519  {
520  EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
521  call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
522  }
523 
526  template<typename OtherDerived>
527  EIGEN_DEVICE_FUNC
528  void swap(MatrixBase<OtherDerived> const & other)
529  {
530  EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
531  call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
532  }
533 
534  template<typename RhsType, typename DstType>
535  EIGEN_DEVICE_FUNC
536  EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const {
537  if(!internal::is_same_dense(dst,rhs))
538  dst = rhs;
539  this->solveInPlace(dst);
540  }
541 
542  template<typename ProductType>
543  EIGEN_DEVICE_FUNC
544  EIGEN_STRONG_INLINE TriangularViewType& _assignProduct(const ProductType& prod, const Scalar& alpha, bool beta);
545  protected:
546  EIGEN_DEFAULT_COPY_CONSTRUCTOR(TriangularViewImpl)
547  EIGEN_DEFAULT_EMPTY_CONSTRUCTOR_AND_DESTRUCTOR(TriangularViewImpl)
548 
549 };
550 
551 /***************************************************************************
552 * Implementation of triangular evaluation/assignment
553 ***************************************************************************/
554 
555 #ifndef EIGEN_PARSED_BY_DOXYGEN
556 // FIXME should we keep that possibility
557 template<typename MatrixType, unsigned int Mode>
558 template<typename OtherDerived>
559 inline TriangularView<MatrixType, Mode>&
560 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const MatrixBase<OtherDerived>& other)
561 {
562  internal::call_assignment_no_alias(derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
563  return derived();
564 }
565 
566 // FIXME should we keep that possibility
567 template<typename MatrixType, unsigned int Mode>
568 template<typename OtherDerived>
569 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const MatrixBase<OtherDerived>& other)
570 {
571  internal::call_assignment_no_alias(derived(), other.template triangularView<Mode>());
572 }
573 
574 
575 
576 template<typename MatrixType, unsigned int Mode>
577 template<typename OtherDerived>
578 inline TriangularView<MatrixType, Mode>&
579 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const TriangularBase<OtherDerived>& other)
580 {
581  eigen_assert(Mode == int(OtherDerived::Mode));
582  internal::call_assignment(derived(), other.derived());
583  return derived();
584 }
585 
586 template<typename MatrixType, unsigned int Mode>
587 template<typename OtherDerived>
588 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const TriangularBase<OtherDerived>& other)
589 {
590  eigen_assert(Mode == int(OtherDerived::Mode));
591  internal::call_assignment_no_alias(derived(), other.derived());
592 }
593 #endif
594 
595 /***************************************************************************
596 * Implementation of TriangularBase methods
597 ***************************************************************************/
598 
601 template<typename Derived>
602 template<typename DenseDerived>
604 {
605  evalToLazy(other.derived());
606 }
607 
608 /***************************************************************************
609 * Implementation of TriangularView methods
610 ***************************************************************************/
611 
612 /***************************************************************************
613 * Implementation of MatrixBase methods
614 ***************************************************************************/
615 
627 template<typename Derived>
628 template<unsigned int Mode>
629 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
631 {
632  return typename TriangularViewReturnType<Mode>::Type(derived());
633 }
634 
636 template<typename Derived>
637 template<unsigned int Mode>
638 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
640 {
641  return typename ConstTriangularViewReturnType<Mode>::Type(derived());
642 }
643 
649 template<typename Derived>
650 bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const
651 {
652  RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
653  for(Index j = 0; j < cols(); ++j)
654  {
655  Index maxi = numext::mini(j, rows()-1);
656  for(Index i = 0; i <= maxi; ++i)
657  {
658  RealScalar absValue = numext::abs(coeff(i,j));
659  if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
660  }
661  }
662  RealScalar threshold = maxAbsOnUpperPart * prec;
663  for(Index j = 0; j < cols(); ++j)
664  for(Index i = j+1; i < rows(); ++i)
665  if(numext::abs(coeff(i, j)) > threshold) return false;
666  return true;
667 }
668 
674 template<typename Derived>
675 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const
676 {
677  RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
678  for(Index j = 0; j < cols(); ++j)
679  for(Index i = j; i < rows(); ++i)
680  {
681  RealScalar absValue = numext::abs(coeff(i,j));
682  if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
683  }
684  RealScalar threshold = maxAbsOnLowerPart * prec;
685  for(Index j = 1; j < cols(); ++j)
686  {
687  Index maxi = numext::mini(j, rows()-1);
688  for(Index i = 0; i < maxi; ++i)
689  if(numext::abs(coeff(i, j)) > threshold) return false;
690  }
691  return true;
692 }
693 
694 
695 /***************************************************************************
696 ****************************************************************************
697 * Evaluators and Assignment of triangular expressions
698 ***************************************************************************
699 ***************************************************************************/
700 
701 namespace internal {
702 
703 
704 // TODO currently a triangular expression has the form TriangularView<.,.>
705 // in the future triangular-ness should be defined by the expression traits
706 // such that Transpose<TriangularView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
707 template<typename MatrixType, unsigned int Mode>
708 struct evaluator_traits<TriangularView<MatrixType,Mode> >
709 {
710  typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
711  typedef typename glue_shapes<typename evaluator_traits<MatrixType>::Shape, TriangularShape>::type Shape;
712 };
713 
714 template<typename MatrixType, unsigned int Mode>
715 struct unary_evaluator<TriangularView<MatrixType,Mode>, IndexBased>
716  : evaluator<typename internal::remove_all<MatrixType>::type>
717 {
718  typedef TriangularView<MatrixType,Mode> XprType;
719  typedef evaluator<typename internal::remove_all<MatrixType>::type> Base;
720  unary_evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {}
721 };
722 
723 // Additional assignment kinds:
724 struct Triangular2Triangular {};
725 struct Triangular2Dense {};
726 struct Dense2Triangular {};
727 
728 
729 template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop;
730 
731 
737 template<int UpLo, int Mode, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version = Specialized>
738 class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version>
739 {
740 protected:
741  typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base;
742  typedef typename Base::DstXprType DstXprType;
743  typedef typename Base::SrcXprType SrcXprType;
744  using Base::m_dst;
745  using Base::m_src;
746  using Base::m_functor;
747 public:
748 
749  typedef typename Base::DstEvaluatorType DstEvaluatorType;
750  typedef typename Base::SrcEvaluatorType SrcEvaluatorType;
751  typedef typename Base::Scalar Scalar;
752  typedef typename Base::AssignmentTraits AssignmentTraits;
753 
754 
755  EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr)
756  : Base(dst, src, func, dstExpr)
757  {}
758 
759 #ifdef EIGEN_INTERNAL_DEBUGGING
760  EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col)
761  {
762  eigen_internal_assert(row!=col);
763  Base::assignCoeff(row,col);
764  }
765 #else
766  using Base::assignCoeff;
767 #endif
768 
769  EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id)
770  {
771  if(Mode==UnitDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(1));
772  else if(Mode==ZeroDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(0));
773  else if(Mode==0) Base::assignCoeff(id,id);
774  }
775 
776  EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index row, Index col)
777  {
778  eigen_internal_assert(row!=col);
779  if(SetOpposite)
780  m_functor.assignCoeff(m_dst.coeffRef(row,col), Scalar(0));
781  }
782 };
783 
784 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType, typename Functor>
785 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
786 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src, const Functor &func)
787 {
788  typedef evaluator<DstXprType> DstEvaluatorType;
789  typedef evaluator<SrcXprType> SrcEvaluatorType;
790 
791  SrcEvaluatorType srcEvaluator(src);
792 
793  Index dstRows = src.rows();
794  Index dstCols = src.cols();
795  if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
796  dst.resize(dstRows, dstCols);
797  DstEvaluatorType dstEvaluator(dst);
798 
799  typedef triangular_dense_assignment_kernel< Mode&(Lower|Upper),Mode&(UnitDiag|ZeroDiag|SelfAdjoint),SetOpposite,
800  DstEvaluatorType,SrcEvaluatorType,Functor> Kernel;
801  Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived());
802 
803  enum {
804  unroll = DstXprType::SizeAtCompileTime != Dynamic
805  && SrcEvaluatorType::CoeffReadCost < HugeCost
806  && DstXprType::SizeAtCompileTime * (DstEvaluatorType::CoeffReadCost+SrcEvaluatorType::CoeffReadCost) / 2 <= EIGEN_UNROLLING_LIMIT
807  };
808 
809  triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, SetOpposite>::run(kernel);
810 }
811 
812 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType>
813 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
814 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src)
815 {
816  call_triangular_assignment_loop<Mode,SetOpposite>(dst, src, internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>());
817 }
818 
819 template<> struct AssignmentKind<TriangularShape,TriangularShape> { typedef Triangular2Triangular Kind; };
820 template<> struct AssignmentKind<DenseShape,TriangularShape> { typedef Triangular2Dense Kind; };
821 template<> struct AssignmentKind<TriangularShape,DenseShape> { typedef Dense2Triangular Kind; };
822 
823 
824 template< typename DstXprType, typename SrcXprType, typename Functor>
825 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Triangular>
826 {
827  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
828  {
829  eigen_assert(int(DstXprType::Mode) == int(SrcXprType::Mode));
830 
831  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
832  }
833 };
834 
835 template< typename DstXprType, typename SrcXprType, typename Functor>
836 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Dense>
837 {
838  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
839  {
840  call_triangular_assignment_loop<SrcXprType::Mode, (SrcXprType::Mode&SelfAdjoint)==0>(dst, src, func);
841  }
842 };
843 
844 template< typename DstXprType, typename SrcXprType, typename Functor>
845 struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular>
846 {
847  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
848  {
849  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
850  }
851 };
852 
853 
854 template<typename Kernel, unsigned int Mode, int UnrollCount, bool SetOpposite>
855 struct triangular_assignment_loop
856 {
857  // FIXME: this is not very clean, perhaps this information should be provided by the kernel?
858  typedef typename Kernel::DstEvaluatorType DstEvaluatorType;
859  typedef typename DstEvaluatorType::XprType DstXprType;
860 
861  enum {
862  col = (UnrollCount-1) / DstXprType::RowsAtCompileTime,
863  row = (UnrollCount-1) % DstXprType::RowsAtCompileTime
864  };
865 
866  typedef typename Kernel::Scalar Scalar;
867 
868  EIGEN_DEVICE_FUNC
869  static inline void run(Kernel &kernel)
870  {
871  triangular_assignment_loop<Kernel, Mode, UnrollCount-1, SetOpposite>::run(kernel);
872 
873  if(row==col)
874  kernel.assignDiagonalCoeff(row);
875  else if( ((Mode&Lower) && row>col) || ((Mode&Upper) && row<col) )
876  kernel.assignCoeff(row,col);
877  else if(SetOpposite)
878  kernel.assignOppositeCoeff(row,col);
879  }
880 };
881 
882 // prevent buggy user code from causing an infinite recursion
883 template<typename Kernel, unsigned int Mode, bool SetOpposite>
884 struct triangular_assignment_loop<Kernel, Mode, 0, SetOpposite>
885 {
886  EIGEN_DEVICE_FUNC
887  static inline void run(Kernel &) {}
888 };
889 
890 
891 
892 // TODO: experiment with a recursive assignment procedure splitting the current
893 // triangular part into one rectangular and two triangular parts.
894 
895 
896 template<typename Kernel, unsigned int Mode, bool SetOpposite>
897 struct triangular_assignment_loop<Kernel, Mode, Dynamic, SetOpposite>
898 {
899  typedef typename Kernel::Scalar Scalar;
900  EIGEN_DEVICE_FUNC
901  static inline void run(Kernel &kernel)
902  {
903  for(Index j = 0; j < kernel.cols(); ++j)
904  {
905  Index maxi = numext::mini(j, kernel.rows());
906  Index i = 0;
907  if (((Mode&Lower) && SetOpposite) || (Mode&Upper))
908  {
909  for(; i < maxi; ++i)
910  if(Mode&Upper) kernel.assignCoeff(i, j);
911  else kernel.assignOppositeCoeff(i, j);
912  }
913  else
914  i = maxi;
915 
916  if(i<kernel.rows()) // then i==j
917  kernel.assignDiagonalCoeff(i++);
918 
919  if (((Mode&Upper) && SetOpposite) || (Mode&Lower))
920  {
921  for(; i < kernel.rows(); ++i)
922  if(Mode&Lower) kernel.assignCoeff(i, j);
923  else kernel.assignOppositeCoeff(i, j);
924  }
925  }
926  }
927 };
928 
929 } // end namespace internal
930 
933 template<typename Derived>
934 template<typename DenseDerived>
936 {
937  other.derived().resize(this->rows(), this->cols());
938  internal::call_triangular_assignment_loop<Derived::Mode,(Derived::Mode&SelfAdjoint)==0 /* SetOpposite */>(other.derived(), derived().nestedExpression());
939 }
940 
941 namespace internal {
942 
943 // Triangular = Product
944 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
945 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
946 {
947  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
948  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,typename SrcXprType::Scalar> &)
949  {
950  Index dstRows = src.rows();
951  Index dstCols = src.cols();
952  if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
953  dst.resize(dstRows, dstCols);
954 
955  dst._assignProduct(src, 1, 0);
956  }
957 };
958 
959 // Triangular += Product
960 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
961 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::add_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::add_assign_op<Scalar,typename SrcXprType::Scalar> &)
965  {
966  dst._assignProduct(src, 1, 1);
967  }
968 };
969 
970 // Triangular -= Product
971 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
972 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::sub_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
973 {
974  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
975  static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,typename SrcXprType::Scalar> &)
976  {
977  dst._assignProduct(src, -1, 1);
978  }
979 };
980 
981 } // end namespace internal
982 
983 } // end namespace Eigen
984 
985 #endif // EIGEN_TRIANGULARMATRIX_H
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::operator*
friend const Product< OtherDerived, TriangularViewType > operator*(const MatrixBase< OtherDerived > &lhs, const TriangularViewImpl &rhs)
Definition: TriangularMatrix.h:460
Eigen::HugeCost
const int HugeCost
Definition: Constants.h:39
Eigen
Namespace containing all symbols from the Eigen library.
Definition: Core:309
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::swap
void swap(TriangularBase< OtherDerived > &other)
Definition: TriangularMatrix.h:515
Eigen::TriangularView::transpose
TransposeReturnType transpose()
Definition: TriangularMatrix.h:252
Eigen::EigenBase::Index
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:38
Eigen::TriangularView::adjoint
const AdjointReturnType adjoint() const
Definition: TriangularMatrix.h:246
Eigen::TriangularView::selfadjointView
const SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView() const
Definition: TriangularMatrix.h:297
Eigen::EigenBase
Definition: EigenBase.h:30
Eigen::TriangularView::conjugate
const ConjugateReturnType conjugate() const
Definition: TriangularMatrix.h:240
Eigen::DenseCoeffsBase< Derived, DirectWriteAccessors >::derived
Derived & derived()
Definition: EigenBase.h:45
Eigen::UnitUpper
@ UnitUpper
Definition: Constants.h:214
Eigen::Upper
@ Upper
Definition: Constants.h:206
Eigen::TriangularView::determinant
Scalar determinant() const
Definition: TriangularMatrix.h:307
Eigen::TriangularView::nestedExpression
NestedExpression & nestedExpression()
Definition: TriangularMatrix.h:235
Eigen::SelfAdjointView
Expression of a selfadjoint matrix from a triangular part of a dense matrix.
Definition: SelfAdjointView.h:51
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::coeffRef
Scalar & coeffRef(Index row, Index col)
Definition: TriangularMatrix.h:414
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::operator*
const Product< TriangularViewType, OtherDerived > operator*(const MatrixBase< OtherDerived > &rhs) const
Definition: TriangularMatrix.h:451
Eigen::DirectAccessBit
const unsigned int DirectAccessBit
Definition: Constants.h:150
Eigen::PacketAccessBit
const unsigned int PacketAccessBit
Definition: Constants.h:89
Eigen::StrictlyUpper
@ StrictlyUpper
Definition: Constants.h:218
Eigen::TriangularBase::evalTo
void evalTo(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:603
Eigen::TriangularBase::copyCoeff
void copyCoeff(Index row, Index col, Other &other)
Definition: TriangularMatrix.h:84
Eigen::LvalueBit
const unsigned int LvalueBit
Definition: Constants.h:139
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::setOnes
TriangularViewType & setOnes()
Definition: TriangularMatrix.h:398
Eigen::TriangularBase
Base class for triangular part in a matrix.
Definition: TriangularMatrix.h:28
Eigen::ZeroDiag
@ ZeroDiag
Definition: Constants.h:210
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:21
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::fill
void fill(const Scalar &value)
Definition: TriangularMatrix.h:388
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::operator+=
TriangularViewType & operator+=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:367
Eigen::Product
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:75
Eigen::StrictlyLower
@ StrictlyLower
Definition: Constants.h:216
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::setZero
TriangularViewType & setZero()
Definition: TriangularMatrix.h:395
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::setConstant
TriangularViewType & setConstant(const Scalar &value)
Definition: TriangularMatrix.h:391
Eigen::Lower
@ Lower
Definition: Constants.h:204
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::coeff
Scalar coeff(Index row, Index col) const
Definition: TriangularMatrix.h:404
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::outerStride
Index outerStride() const
Definition: TriangularMatrix.h:358
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::operator-=
TriangularViewType & operator-=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:374
Eigen::TriangularView::cols
Index cols() const
Definition: TriangularMatrix.h:227
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::solve
const internal::triangular_solve_retval< Side, TriangularViewType, Other > solve(const MatrixBase< Other > &other) const
Eigen::TriangularView::rows
Index rows() const
Definition: TriangularMatrix.h:224
Eigen::EigenBase::derived
Derived & derived()
Definition: EigenBase.h:45
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::operator=
TriangularViewType & operator=(const TriangularBase< OtherDerived > &other)
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::swap
void swap(MatrixBase< OtherDerived > const &other)
Definition: TriangularMatrix.h:528
Eigen::MatrixBase::isLowerTriangular
bool isLowerTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:675
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::innerStride
Index innerStride() const
Definition: TriangularMatrix.h:362
Eigen::LinearAccessBit
const unsigned int LinearAccessBit
Definition: Constants.h:125
Eigen::Solve
Pseudo expression representing a solving operation.
Definition: Solve.h:63
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::solveInPlace
void solveInPlace(const MatrixBase< OtherDerived > &other) const
Eigen::DenseBase
Base class for all dense matrices, vectors, and arrays.
Definition: DenseBase.h:47
Eigen::TriangularView::nestedExpression
const NestedExpression & nestedExpression() const
Definition: TriangularMatrix.h:231
Eigen::TriangularView::transpose
const ConstTransposeReturnType transpose() const
Definition: TriangularMatrix.h:262
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::operator/=
TriangularViewType & operator/=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:384
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::operator=
TriangularViewType & operator=(const MatrixBase< OtherDerived > &other)
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
Eigen::TriangularView::selfadjointView
SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView()
Definition: TriangularMatrix.h:289
Eigen::TriangularViewImpl< _MatrixType, _Mode, Dense >::operator*=
TriangularViewType & operator*=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:381
Eigen::TriangularView
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:188
Eigen::TriangularBase::evalToLazy
void evalToLazy(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:935
Eigen::UnitLower
@ UnitLower
Definition: Constants.h:212
Eigen::SelfAdjoint
@ SelfAdjoint
Definition: Constants.h:220
Eigen::MatrixBase::isUpperTriangular
bool isUpperTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:650
Eigen::UnitDiag
@ UnitDiag
Definition: Constants.h:208
Eigen::TriangularBase::SizeAtCompileTime
@ SizeAtCompileTime
Definition: TriangularMatrix.h:38
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Eigen::Dense
Definition: Constants.h:491