Eigen  3.3.0
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
114  void evalToLazy(MatrixBase<DenseDerived> &other) const;
115 
116  EIGEN_DEVICE_FUNC
117  DenseMatrixType toDenseMatrix() const
118  {
119  DenseMatrixType res(rows(), cols());
120  evalToLazy(res);
121  return res;
122  }
123 
124  protected:
125 
126  void check_coordinates(Index row, Index col) const
127  {
128  EIGEN_ONLY_USED_FOR_DEBUG(row);
129  EIGEN_ONLY_USED_FOR_DEBUG(col);
130  eigen_assert(col>=0 && col<cols() && row>=0 && row<rows());
131  const int mode = int(Mode) & ~SelfAdjoint;
132  EIGEN_ONLY_USED_FOR_DEBUG(mode);
133  eigen_assert((mode==Upper && col>=row)
134  || (mode==Lower && col<=row)
135  || ((mode==StrictlyUpper || mode==UnitUpper) && col>row)
136  || ((mode==StrictlyLower || mode==UnitLower) && col<row));
137  }
138 
139  #ifdef EIGEN_INTERNAL_DEBUGGING
140  void check_coordinates_internal(Index row, Index col) const
141  {
142  check_coordinates(row, col);
143  }
144  #else
145  void check_coordinates_internal(Index , Index ) const {}
146  #endif
147 
148 };
149 
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  using Base::operator=;
221  TriangularView& operator=(const TriangularView &other)
222  { return Base::operator=(other); }
223 
225  EIGEN_DEVICE_FUNC
226  inline Index rows() const { return m_matrix.rows(); }
228  EIGEN_DEVICE_FUNC
229  inline Index cols() const { 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 
247  EIGEN_DEVICE_FUNC
248  inline const AdjointReturnType adjoint() const
249  { return AdjointReturnType(m_matrix.adjoint()); }
250 
253  EIGEN_DEVICE_FUNC
254  inline TransposeReturnType transpose()
255  {
256  EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
257  typename MatrixType::TransposeReturnType tmp(m_matrix);
258  return TransposeReturnType(tmp);
259  }
260 
263  EIGEN_DEVICE_FUNC
264  inline const ConstTransposeReturnType transpose() const
265  {
266  return ConstTransposeReturnType(m_matrix.transpose());
267  }
268 
269  template<typename Other>
270  EIGEN_DEVICE_FUNC
271  inline const Solve<TriangularView, Other>
272  solve(const MatrixBase<Other>& other) const
273  { return Solve<TriangularView, Other>(*this, other.derived()); }
274 
275  // workaround MSVC ICE
276  #if EIGEN_COMP_MSVC
277  template<int Side, typename Other>
278  EIGEN_DEVICE_FUNC
279  inline const internal::triangular_solve_retval<Side,TriangularView, Other>
280  solve(const MatrixBase<Other>& other) const
281  { return Base::template solve<Side>(other); }
282  #else
283  using Base::solve;
284  #endif
285 
290  EIGEN_DEVICE_FUNC
292  {
293  EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
295  }
296 
298  EIGEN_DEVICE_FUNC
300  {
301  EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
303  }
304 
305 
308  EIGEN_DEVICE_FUNC
309  Scalar determinant() const
310  {
311  if (Mode & UnitDiag)
312  return 1;
313  else if (Mode & ZeroDiag)
314  return 0;
315  else
316  return m_matrix.diagonal().prod();
317  }
318 
319  protected:
320 
321  MatrixTypeNested m_matrix;
322 };
323 
333 template<typename _MatrixType, unsigned int _Mode> class TriangularViewImpl<_MatrixType,_Mode,Dense>
334  : public TriangularBase<TriangularView<_MatrixType, _Mode> >
335 {
336  public:
337 
340  typedef typename internal::traits<TriangularViewType>::Scalar Scalar;
341 
342  typedef _MatrixType MatrixType;
343  typedef typename MatrixType::PlainObject DenseMatrixType;
344  typedef DenseMatrixType PlainObject;
345 
346  public:
347  using Base::evalToLazy;
348  using Base::derived;
349 
350  typedef typename internal::traits<TriangularViewType>::StorageKind StorageKind;
351 
352  enum {
353  Mode = _Mode,
354  Flags = internal::traits<TriangularViewType>::Flags
355  };
356 
359  EIGEN_DEVICE_FUNC
360  inline Index outerStride() const { return derived().nestedExpression().outerStride(); }
363  EIGEN_DEVICE_FUNC
364  inline Index innerStride() const { return derived().nestedExpression().innerStride(); }
365 
367  template<typename Other>
368  EIGEN_DEVICE_FUNC
369  TriangularViewType& operator+=(const DenseBase<Other>& other) {
370  internal::call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op<Scalar,typename Other::Scalar>());
371  return derived();
372  }
374  template<typename Other>
375  EIGEN_DEVICE_FUNC
376  TriangularViewType& operator-=(const DenseBase<Other>& other) {
377  internal::call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op<Scalar,typename Other::Scalar>());
378  return derived();
379  }
380 
382  EIGEN_DEVICE_FUNC
383  TriangularViewType& operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() * other; }
385  EIGEN_DEVICE_FUNC
386  TriangularViewType& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() / other; }
387 
389  EIGEN_DEVICE_FUNC
390  void fill(const Scalar& value) { setConstant(value); }
392  EIGEN_DEVICE_FUNC
393  TriangularViewType& setConstant(const Scalar& value)
394  { return *this = MatrixType::Constant(derived().rows(), derived().cols(), value); }
396  EIGEN_DEVICE_FUNC
397  TriangularViewType& setZero() { return setConstant(Scalar(0)); }
399  EIGEN_DEVICE_FUNC
400  TriangularViewType& setOnes() { return setConstant(Scalar(1)); }
401 
405  EIGEN_DEVICE_FUNC
406  inline Scalar coeff(Index row, Index col) const
407  {
408  Base::check_coordinates_internal(row, col);
409  return derived().nestedExpression().coeff(row, col);
410  }
411 
415  EIGEN_DEVICE_FUNC
416  inline Scalar& coeffRef(Index row, Index col)
417  {
418  EIGEN_STATIC_ASSERT_LVALUE(TriangularViewType);
419  Base::check_coordinates_internal(row, col);
420  return derived().nestedExpression().coeffRef(row, col);
421  }
422 
424  template<typename OtherDerived>
425  EIGEN_DEVICE_FUNC
426  TriangularViewType& operator=(const TriangularBase<OtherDerived>& other);
427 
429  template<typename OtherDerived>
430  EIGEN_DEVICE_FUNC
431  TriangularViewType& operator=(const MatrixBase<OtherDerived>& other);
432 
433 #ifndef EIGEN_PARSED_BY_DOXYGEN
434  EIGEN_DEVICE_FUNC
435  TriangularViewType& operator=(const TriangularViewImpl& other)
436  { return *this = other.derived().nestedExpression(); }
437 
439  template<typename OtherDerived>
440  EIGEN_DEVICE_FUNC
441  void lazyAssign(const TriangularBase<OtherDerived>& other);
442 
444  template<typename OtherDerived>
445  EIGEN_DEVICE_FUNC
446  void lazyAssign(const MatrixBase<OtherDerived>& other);
447 #endif
448 
450  template<typename OtherDerived>
451  EIGEN_DEVICE_FUNC
454  {
455  return Product<TriangularViewType,OtherDerived>(derived(), rhs.derived());
456  }
457 
459  template<typename OtherDerived> friend
460  EIGEN_DEVICE_FUNC
462  operator*(const MatrixBase<OtherDerived>& lhs, const TriangularViewImpl& rhs)
463  {
464  return Product<OtherDerived,TriangularViewType>(lhs.derived(),rhs.derived());
465  }
466 
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 
500  template<int Side, typename OtherDerived>
501  EIGEN_DEVICE_FUNC
502  void solveInPlace(const MatrixBase<OtherDerived>& other) const;
503 
504  template<typename OtherDerived>
505  EIGEN_DEVICE_FUNC
506  void solveInPlace(const MatrixBase<OtherDerived>& other) const
507  { return solveInPlace<OnTheLeft>(other); }
508 
510  template<typename OtherDerived>
511  EIGEN_DEVICE_FUNC
512 #ifdef EIGEN_PARSED_BY_DOXYGEN
514 #else
515  void swap(TriangularBase<OtherDerived> const & other)
516 #endif
517  {
518  EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
519  call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
520  }
521 
524  template<typename OtherDerived>
525  EIGEN_DEVICE_FUNC
526  void swap(MatrixBase<OtherDerived> const & other)
527  {
528  EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
529  call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
530  }
531 
532  template<typename RhsType, typename DstType>
533  EIGEN_DEVICE_FUNC
534  EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const {
535  if(!internal::is_same_dense(dst,rhs))
536  dst = rhs;
537  this->solveInPlace(dst);
538  }
539 
540  template<typename ProductType>
541  EIGEN_DEVICE_FUNC
542  EIGEN_STRONG_INLINE TriangularViewType& _assignProduct(const ProductType& prod, const Scalar& alpha);
543 };
544 
545 /***************************************************************************
546 * Implementation of triangular evaluation/assignment
547 ***************************************************************************/
548 
549 // FIXME should we keep that possibility
550 template<typename MatrixType, unsigned int Mode>
551 template<typename OtherDerived>
553 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const MatrixBase<OtherDerived>& other)
554 {
555  internal::call_assignment_no_alias(derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
556  return derived();
557 }
558 
559 // FIXME should we keep that possibility
560 template<typename MatrixType, unsigned int Mode>
561 template<typename OtherDerived>
562 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const MatrixBase<OtherDerived>& other)
563 {
564  internal::call_assignment_no_alias(derived(), other.template triangularView<Mode>());
565 }
566 
567 
568 
569 template<typename MatrixType, unsigned int Mode>
570 template<typename OtherDerived>
572 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const TriangularBase<OtherDerived>& other)
573 {
574  eigen_assert(Mode == int(OtherDerived::Mode));
575  internal::call_assignment(derived(), other.derived());
576  return derived();
577 }
578 
579 template<typename MatrixType, unsigned int Mode>
580 template<typename OtherDerived>
581 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const TriangularBase<OtherDerived>& other)
582 {
583  eigen_assert(Mode == int(OtherDerived::Mode));
584  internal::call_assignment_no_alias(derived(), other.derived());
585 }
586 
587 /***************************************************************************
588 * Implementation of TriangularBase methods
589 ***************************************************************************/
590 
593 template<typename Derived>
594 template<typename DenseDerived>
596 {
597  evalToLazy(other.derived());
598 }
599 
600 /***************************************************************************
601 * Implementation of TriangularView methods
602 ***************************************************************************/
603 
604 /***************************************************************************
605 * Implementation of MatrixBase methods
606 ***************************************************************************/
607 
619 template<typename Derived>
620 template<unsigned int Mode>
621 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
623 {
624  return typename TriangularViewReturnType<Mode>::Type(derived());
625 }
626 
628 template<typename Derived>
629 template<unsigned int Mode>
630 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
632 {
633  return typename ConstTriangularViewReturnType<Mode>::Type(derived());
634 }
635 
641 template<typename Derived>
642 bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const
643 {
644  RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
645  for(Index j = 0; j < cols(); ++j)
646  {
647  Index maxi = numext::mini(j, rows()-1);
648  for(Index i = 0; i <= maxi; ++i)
649  {
650  RealScalar absValue = numext::abs(coeff(i,j));
651  if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
652  }
653  }
654  RealScalar threshold = maxAbsOnUpperPart * prec;
655  for(Index j = 0; j < cols(); ++j)
656  for(Index i = j+1; i < rows(); ++i)
657  if(numext::abs(coeff(i, j)) > threshold) return false;
658  return true;
659 }
660 
666 template<typename Derived>
667 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const
668 {
669  RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
670  for(Index j = 0; j < cols(); ++j)
671  for(Index i = j; i < rows(); ++i)
672  {
673  RealScalar absValue = numext::abs(coeff(i,j));
674  if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
675  }
676  RealScalar threshold = maxAbsOnLowerPart * prec;
677  for(Index j = 1; j < cols(); ++j)
678  {
679  Index maxi = numext::mini(j, rows()-1);
680  for(Index i = 0; i < maxi; ++i)
681  if(numext::abs(coeff(i, j)) > threshold) return false;
682  }
683  return true;
684 }
685 
686 
687 /***************************************************************************
688 ****************************************************************************
689 * Evaluators and Assignment of triangular expressions
690 ***************************************************************************
691 ***************************************************************************/
692 
693 namespace internal {
694 
695 
696 // TODO currently a triangular expression has the form TriangularView<.,.>
697 // in the future triangular-ness should be defined by the expression traits
698 // such that Transpose<TriangularView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
699 template<typename MatrixType, unsigned int Mode>
700 struct evaluator_traits<TriangularView<MatrixType,Mode> >
701 {
702  typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
703  typedef typename glue_shapes<typename evaluator_traits<MatrixType>::Shape, TriangularShape>::type Shape;
704 };
705 
706 template<typename MatrixType, unsigned int Mode>
707 struct unary_evaluator<TriangularView<MatrixType,Mode>, IndexBased>
708  : evaluator<typename internal::remove_all<MatrixType>::type>
709 {
710  typedef TriangularView<MatrixType,Mode> XprType;
711  typedef evaluator<typename internal::remove_all<MatrixType>::type> Base;
712  unary_evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {}
713 };
714 
715 // Additional assignment kinds:
716 struct Triangular2Triangular {};
717 struct Triangular2Dense {};
718 struct Dense2Triangular {};
719 
720 
721 template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop;
722 
723 
729 template<int UpLo, int Mode, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version = Specialized>
730 class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version>
731 {
732 protected:
733  typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base;
734  typedef typename Base::DstXprType DstXprType;
735  typedef typename Base::SrcXprType SrcXprType;
736  using Base::m_dst;
737  using Base::m_src;
738  using Base::m_functor;
739 public:
740 
741  typedef typename Base::DstEvaluatorType DstEvaluatorType;
742  typedef typename Base::SrcEvaluatorType SrcEvaluatorType;
743  typedef typename Base::Scalar Scalar;
744  typedef typename Base::AssignmentTraits AssignmentTraits;
745 
746 
747  EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr)
748  : Base(dst, src, func, dstExpr)
749  {}
750 
751 #ifdef EIGEN_INTERNAL_DEBUGGING
752  EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col)
753  {
754  eigen_internal_assert(row!=col);
755  Base::assignCoeff(row,col);
756  }
757 #else
758  using Base::assignCoeff;
759 #endif
760 
761  EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id)
762  {
763  if(Mode==UnitDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(1));
764  else if(Mode==ZeroDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(0));
765  else if(Mode==0) Base::assignCoeff(id,id);
766  }
767 
768  EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index row, Index col)
769  {
770  eigen_internal_assert(row!=col);
771  if(SetOpposite)
772  m_functor.assignCoeff(m_dst.coeffRef(row,col), Scalar(0));
773  }
774 };
775 
776 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType, typename Functor>
777 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
778 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src, const Functor &func)
779 {
780  typedef evaluator<DstXprType> DstEvaluatorType;
781  typedef evaluator<SrcXprType> SrcEvaluatorType;
782 
783  SrcEvaluatorType srcEvaluator(src);
784 
785  Index dstRows = src.rows();
786  Index dstCols = src.cols();
787  if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
788  dst.resize(dstRows, dstCols);
789  DstEvaluatorType dstEvaluator(dst);
790 
791  typedef triangular_dense_assignment_kernel< Mode&(Lower|Upper),Mode&(UnitDiag|ZeroDiag|SelfAdjoint),SetOpposite,
792  DstEvaluatorType,SrcEvaluatorType,Functor> Kernel;
793  Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived());
794 
795  enum {
796  unroll = DstXprType::SizeAtCompileTime != Dynamic
797  && SrcEvaluatorType::CoeffReadCost < HugeCost
798  && DstXprType::SizeAtCompileTime * (DstEvaluatorType::CoeffReadCost+SrcEvaluatorType::CoeffReadCost) / 2 <= EIGEN_UNROLLING_LIMIT
799  };
800 
801  triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, SetOpposite>::run(kernel);
802 }
803 
804 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType>
805 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
806 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src)
807 {
808  call_triangular_assignment_loop<Mode,SetOpposite>(dst, src, internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>());
809 }
810 
811 template<> struct AssignmentKind<TriangularShape,TriangularShape> { typedef Triangular2Triangular Kind; };
812 template<> struct AssignmentKind<DenseShape,TriangularShape> { typedef Triangular2Dense Kind; };
813 template<> struct AssignmentKind<TriangularShape,DenseShape> { typedef Dense2Triangular Kind; };
814 
815 
816 template< typename DstXprType, typename SrcXprType, typename Functor>
817 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Triangular>
818 {
819  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
820  {
821  eigen_assert(int(DstXprType::Mode) == int(SrcXprType::Mode));
822 
823  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
824  }
825 };
826 
827 template< typename DstXprType, typename SrcXprType, typename Functor>
828 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Dense>
829 {
830  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
831  {
832  call_triangular_assignment_loop<SrcXprType::Mode, (SrcXprType::Mode&SelfAdjoint)==0>(dst, src, func);
833  }
834 };
835 
836 template< typename DstXprType, typename SrcXprType, typename Functor>
837 struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular>
838 {
839  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
840  {
841  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
842  }
843 };
844 
845 
846 template<typename Kernel, unsigned int Mode, int UnrollCount, bool SetOpposite>
847 struct triangular_assignment_loop
848 {
849  // FIXME: this is not very clean, perhaps this information should be provided by the kernel?
850  typedef typename Kernel::DstEvaluatorType DstEvaluatorType;
851  typedef typename DstEvaluatorType::XprType DstXprType;
852 
853  enum {
854  col = (UnrollCount-1) / DstXprType::RowsAtCompileTime,
855  row = (UnrollCount-1) % DstXprType::RowsAtCompileTime
856  };
857 
858  typedef typename Kernel::Scalar Scalar;
859 
860  EIGEN_DEVICE_FUNC
861  static inline void run(Kernel &kernel)
862  {
863  triangular_assignment_loop<Kernel, Mode, UnrollCount-1, SetOpposite>::run(kernel);
864 
865  if(row==col)
866  kernel.assignDiagonalCoeff(row);
867  else if( ((Mode&Lower) && row>col) || ((Mode&Upper) && row<col) )
868  kernel.assignCoeff(row,col);
869  else if(SetOpposite)
870  kernel.assignOppositeCoeff(row,col);
871  }
872 };
873 
874 // prevent buggy user code from causing an infinite recursion
875 template<typename Kernel, unsigned int Mode, bool SetOpposite>
876 struct triangular_assignment_loop<Kernel, Mode, 0, SetOpposite>
877 {
878  EIGEN_DEVICE_FUNC
879  static inline void run(Kernel &) {}
880 };
881 
882 
883 
884 // TODO: experiment with a recursive assignment procedure splitting the current
885 // triangular part into one rectangular and two triangular parts.
886 
887 
888 template<typename Kernel, unsigned int Mode, bool SetOpposite>
889 struct triangular_assignment_loop<Kernel, Mode, Dynamic, SetOpposite>
890 {
891  typedef typename Kernel::Scalar Scalar;
892  EIGEN_DEVICE_FUNC
893  static inline void run(Kernel &kernel)
894  {
895  for(Index j = 0; j < kernel.cols(); ++j)
896  {
897  Index maxi = numext::mini(j, kernel.rows());
898  Index i = 0;
899  if (((Mode&Lower) && SetOpposite) || (Mode&Upper))
900  {
901  for(; i < maxi; ++i)
902  if(Mode&Upper) kernel.assignCoeff(i, j);
903  else kernel.assignOppositeCoeff(i, j);
904  }
905  else
906  i = maxi;
907 
908  if(i<kernel.rows()) // then i==j
909  kernel.assignDiagonalCoeff(i++);
910 
911  if (((Mode&Upper) && SetOpposite) || (Mode&Lower))
912  {
913  for(; i < kernel.rows(); ++i)
914  if(Mode&Lower) kernel.assignCoeff(i, j);
915  else kernel.assignOppositeCoeff(i, j);
916  }
917  }
918  }
919 };
920 
921 } // end namespace internal
922 
925 template<typename Derived>
926 template<typename DenseDerived>
928 {
929  other.derived().resize(this->rows(), this->cols());
930  internal::call_triangular_assignment_loop<Derived::Mode,(Derived::Mode&SelfAdjoint)==0 /* SetOpposite */>(other.derived(), derived().nestedExpression());
931 }
932 
933 namespace internal {
934 
935 // Triangular = Product
936 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
937 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
938 {
939  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
940  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,typename SrcXprType::Scalar> &)
941  {
942  Index dstRows = src.rows();
943  Index dstCols = src.cols();
944  if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
945  dst.resize(dstRows, dstCols);
946 
947  dst.setZero();
948  dst._assignProduct(src, 1);
949  }
950 };
951 
952 // Triangular += Product
953 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
954 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::add_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
955 {
956  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
957  static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar,typename SrcXprType::Scalar> &)
958  {
959  dst._assignProduct(src, 1);
960  }
961 };
962 
963 // Triangular -= Product
964 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
965 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::sub_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
966 {
967  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
968  static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,typename SrcXprType::Scalar> &)
969  {
970  dst._assignProduct(src, -1);
971  }
972 };
973 
974 } // end namespace internal
975 
976 } // end namespace Eigen
977 
978 #endif // EIGEN_TRIANGULARMATRIX_H
const Product< TriangularViewType, OtherDerived > operator*(const MatrixBase< OtherDerived > &rhs) const
Definition: TriangularMatrix.h:453
Scalar determinant() const
Definition: TriangularMatrix.h:309
const int HugeCost
Definition: Constants.h:39
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:71
Index outerStride() const
Definition: TriangularMatrix.h:360
NestedExpression & nestedExpression()
Definition: TriangularMatrix.h:237
TriangularViewType & setConstant(const Scalar &value)
Definition: TriangularMatrix.h:393
Definition: Constants.h:216
TransposeReturnType transpose()
Definition: TriangularMatrix.h:254
Base class for triangular part in a matrix.
Definition: TriangularMatrix.h:27
const unsigned int DirectAccessBit
Definition: Constants.h:150
const AdjointReturnType adjoint() const
Definition: TriangularMatrix.h:248
const SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView() const
Definition: TriangularMatrix.h:299
const unsigned int LvalueBit
Definition: Constants.h:139
friend const Product< OtherDerived, TriangularViewType > operator*(const MatrixBase< OtherDerived > &lhs, const TriangularViewImpl &rhs)
Definition: TriangularMatrix.h:462
Namespace containing all symbols from the Eigen library.
Definition: Core:287
TriangularViewType & setOnes()
Definition: TriangularMatrix.h:400
Derived & derived()
Definition: EigenBase.h:44
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:37
TriangularViewType & operator*=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:383
Base class for all dense matrices, vectors, and arrays.
Definition: DenseBase.h:41
const unsigned int PacketAccessBit
Definition: Constants.h:89
const ConjugateReturnType conjugate() const
Definition: TriangularMatrix.h:242
Definition: EigenBase.h:28
bool isUpperTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:642
void copyCoeff(Index row, Index col, Other &other)
Definition: TriangularMatrix.h:84
Definition: Constants.h:204
void swap(MatrixBase< OtherDerived > const &other)
Definition: TriangularMatrix.h:526
Definition: Constants.h:208
Scalar coeff(Index row, Index col) const
Definition: TriangularMatrix.h:406
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
void swap(TriangularBase< OtherDerived > &other)
Definition: TriangularMatrix.h:513
void evalTo(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:595
TriangularViewType & operator/=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:386
Expression of a selfadjoint matrix from a triangular part of a dense matrix.
Definition: SelfAdjointView.h:49
Definition: Constants.h:218
Definition: Constants.h:214
Definition: Constants.h:206
TriangularViewType & setZero()
Definition: TriangularMatrix.h:397
Definition: Constants.h:210
Definition: Constants.h:212
Definition: Eigen_Colamd.h:50
bool isLowerTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:667
Index rows() const
Definition: TriangularMatrix.h:226
Scalar & coeffRef(Index row, Index col)
Definition: TriangularMatrix.h:416
Definition: Constants.h:220
Index cols() const
Definition: TriangularMatrix.h:229
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:186
SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView()
Definition: TriangularMatrix.h:291
TriangularViewType & operator-=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:376
Definition: Constants.h:491
Index innerStride() const
Definition: TriangularMatrix.h:364
const int Dynamic
Definition: Constants.h:21
Pseudo expression representing a solving operation.
Definition: Solve.h:62
const NestedExpression & nestedExpression() const
Definition: TriangularMatrix.h:233
const ConstTransposeReturnType transpose() const
Definition: TriangularMatrix.h:264
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const unsigned int LinearAccessBit
Definition: Constants.h:125
void fill(const Scalar &value)
Definition: TriangularMatrix.h:390
void evalToLazy(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:927
TriangularViewType & operator+=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:369