Eigen  3.2.93
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>
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  using std::abs;
645  RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
646  for(Index j = 0; j < cols(); ++j)
647  {
648  Index maxi = (std::min)(j, rows()-1);
649  for(Index i = 0; i <= maxi; ++i)
650  {
651  RealScalar absValue = abs(coeff(i,j));
652  if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
653  }
654  }
655  RealScalar threshold = maxAbsOnUpperPart * prec;
656  for(Index j = 0; j < cols(); ++j)
657  for(Index i = j+1; i < rows(); ++i)
658  if(abs(coeff(i, j)) > threshold) return false;
659  return true;
660 }
661 
667 template<typename Derived>
668 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const
669 {
670  using std::abs;
671  RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
672  for(Index j = 0; j < cols(); ++j)
673  for(Index i = j; i < rows(); ++i)
674  {
675  RealScalar absValue = abs(coeff(i,j));
676  if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
677  }
678  RealScalar threshold = maxAbsOnLowerPart * prec;
679  for(Index j = 1; j < cols(); ++j)
680  {
681  Index maxi = (std::min)(j, rows()-1);
682  for(Index i = 0; i < maxi; ++i)
683  if(abs(coeff(i, j)) > threshold) return false;
684  }
685  return true;
686 }
687 
688 
689 /***************************************************************************
690 ****************************************************************************
691 * Evaluators and Assignment of triangular expressions
692 ***************************************************************************
693 ***************************************************************************/
694 
695 namespace internal {
696 
697 
698 // TODO currently a triangular expression has the form TriangularView<.,.>
699 // in the future triangular-ness should be defined by the expression traits
700 // such that Transpose<TriangularView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
701 template<typename MatrixType, unsigned int Mode>
702 struct evaluator_traits<TriangularView<MatrixType,Mode> >
703 {
704  typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
705  typedef typename glue_shapes<typename evaluator_traits<MatrixType>::Shape, TriangularShape>::type Shape;
706 };
707 
708 template<typename MatrixType, unsigned int Mode>
709 struct unary_evaluator<TriangularView<MatrixType,Mode>, IndexBased>
710  : evaluator<typename internal::remove_all<MatrixType>::type>
711 {
712  typedef TriangularView<MatrixType,Mode> XprType;
713  typedef evaluator<typename internal::remove_all<MatrixType>::type> Base;
714  unary_evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {}
715 };
716 
717 // Additional assignment kinds:
718 struct Triangular2Triangular {};
719 struct Triangular2Dense {};
720 struct Dense2Triangular {};
721 
722 
723 template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop;
724 
725 
731 template<int UpLo, int Mode, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version = Specialized>
732 class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version>
733 {
734 protected:
735  typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base;
736  typedef typename Base::DstXprType DstXprType;
737  typedef typename Base::SrcXprType SrcXprType;
738  using Base::m_dst;
739  using Base::m_src;
740  using Base::m_functor;
741 public:
742 
743  typedef typename Base::DstEvaluatorType DstEvaluatorType;
744  typedef typename Base::SrcEvaluatorType SrcEvaluatorType;
745  typedef typename Base::Scalar Scalar;
746  typedef typename Base::AssignmentTraits AssignmentTraits;
747 
748 
749  EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr)
750  : Base(dst, src, func, dstExpr)
751  {}
752 
753 #ifdef EIGEN_INTERNAL_DEBUGGING
754  EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col)
755  {
756  eigen_internal_assert(row!=col);
757  Base::assignCoeff(row,col);
758  }
759 #else
760  using Base::assignCoeff;
761 #endif
762 
763  EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id)
764  {
765  if(Mode==UnitDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(1));
766  else if(Mode==ZeroDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(0));
767  else if(Mode==0) Base::assignCoeff(id,id);
768  }
769 
770  EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index row, Index col)
771  {
772  eigen_internal_assert(row!=col);
773  if(SetOpposite)
774  m_functor.assignCoeff(m_dst.coeffRef(row,col), Scalar(0));
775  }
776 };
777 
778 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType, typename Functor>
779 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
780 void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src, const Functor &func)
781 {
782  eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
783 
784  typedef evaluator<DstXprType> DstEvaluatorType;
785  typedef evaluator<SrcXprType> SrcEvaluatorType;
786 
787  DstEvaluatorType dstEvaluator(dst);
788  SrcEvaluatorType srcEvaluator(src);
789 
790  typedef triangular_dense_assignment_kernel< Mode&(Lower|Upper),Mode&(UnitDiag|ZeroDiag|SelfAdjoint),SetOpposite,
791  DstEvaluatorType,SrcEvaluatorType,Functor> Kernel;
792  Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived());
793 
794  enum {
795  unroll = DstXprType::SizeAtCompileTime != Dynamic
796  && SrcEvaluatorType::CoeffReadCost < HugeCost
797  && DstXprType::SizeAtCompileTime * (DstEvaluatorType::CoeffReadCost+SrcEvaluatorType::CoeffReadCost) / 2 <= EIGEN_UNROLLING_LIMIT
798  };
799 
800  triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, SetOpposite>::run(kernel);
801 }
802 
803 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType>
804 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
805 void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src)
806 {
807  call_triangular_assignment_loop<Mode,SetOpposite>(dst, src, internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>());
808 }
809 
810 template<> struct AssignmentKind<TriangularShape,TriangularShape> { typedef Triangular2Triangular Kind; };
811 template<> struct AssignmentKind<DenseShape,TriangularShape> { typedef Triangular2Dense Kind; };
812 template<> struct AssignmentKind<TriangularShape,DenseShape> { typedef Dense2Triangular Kind; };
813 
814 
815 template< typename DstXprType, typename SrcXprType, typename Functor>
816 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Triangular>
817 {
818  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
819  {
820  eigen_assert(int(DstXprType::Mode) == int(SrcXprType::Mode));
821 
822  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
823  }
824 };
825 
826 template< typename DstXprType, typename SrcXprType, typename Functor>
827 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Dense>
828 {
829  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
830  {
831  call_triangular_assignment_loop<SrcXprType::Mode, (SrcXprType::Mode&SelfAdjoint)==0>(dst, src, func);
832  }
833 };
834 
835 template< typename DstXprType, typename SrcXprType, typename Functor>
836 struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular>
837 {
838  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
839  {
840  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
841  }
842 };
843 
844 
845 template<typename Kernel, unsigned int Mode, int UnrollCount, bool SetOpposite>
846 struct triangular_assignment_loop
847 {
848  // FIXME: this is not very clean, perhaps this information should be provided by the kernel?
849  typedef typename Kernel::DstEvaluatorType DstEvaluatorType;
850  typedef typename DstEvaluatorType::XprType DstXprType;
851 
852  enum {
853  col = (UnrollCount-1) / DstXprType::RowsAtCompileTime,
854  row = (UnrollCount-1) % DstXprType::RowsAtCompileTime
855  };
856 
857  typedef typename Kernel::Scalar Scalar;
858 
859  EIGEN_DEVICE_FUNC
860  static inline void run(Kernel &kernel)
861  {
862  triangular_assignment_loop<Kernel, Mode, UnrollCount-1, SetOpposite>::run(kernel);
863 
864  if(row==col)
865  kernel.assignDiagonalCoeff(row);
866  else if( ((Mode&Lower) && row>col) || ((Mode&Upper) && row<col) )
867  kernel.assignCoeff(row,col);
868  else if(SetOpposite)
869  kernel.assignOppositeCoeff(row,col);
870  }
871 };
872 
873 // prevent buggy user code from causing an infinite recursion
874 template<typename Kernel, unsigned int Mode, bool SetOpposite>
875 struct triangular_assignment_loop<Kernel, Mode, 0, SetOpposite>
876 {
877  EIGEN_DEVICE_FUNC
878  static inline void run(Kernel &) {}
879 };
880 
881 
882 
883 // TODO: experiment with a recursive assignment procedure splitting the current
884 // triangular part into one rectangular and two triangular parts.
885 
886 
887 template<typename Kernel, unsigned int Mode, bool SetOpposite>
888 struct triangular_assignment_loop<Kernel, Mode, Dynamic, SetOpposite>
889 {
890  typedef typename Kernel::Scalar Scalar;
891  EIGEN_DEVICE_FUNC
892  static inline void run(Kernel &kernel)
893  {
894  for(Index j = 0; j < kernel.cols(); ++j)
895  {
896  Index maxi = (std::min)(j, kernel.rows());
897  Index i = 0;
898  if (((Mode&Lower) && SetOpposite) || (Mode&Upper))
899  {
900  for(; i < maxi; ++i)
901  if(Mode&Upper) kernel.assignCoeff(i, j);
902  else kernel.assignOppositeCoeff(i, j);
903  }
904  else
905  i = maxi;
906 
907  if(i<kernel.rows()) // then i==j
908  kernel.assignDiagonalCoeff(i++);
909 
910  if (((Mode&Upper) && SetOpposite) || (Mode&Lower))
911  {
912  for(; i < kernel.rows(); ++i)
913  if(Mode&Lower) kernel.assignCoeff(i, j);
914  else kernel.assignOppositeCoeff(i, j);
915  }
916  }
917  }
918 };
919 
920 } // end namespace internal
921 
924 template<typename Derived>
925 template<typename DenseDerived>
927 {
928  other.derived().resize(this->rows(), this->cols());
929  internal::call_triangular_assignment_loop<Derived::Mode,(Derived::Mode&SelfAdjoint)==0 /* SetOpposite */>(other.derived(), derived().nestedExpression());
930 }
931 
932 namespace internal {
933 
934 // Triangular = Product
935 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
936 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
937 {
938  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
939  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,typename SrcXprType::Scalar> &)
940  {
941  dst.setZero();
942  dst._assignProduct(src, 1);
943  }
944 };
945 
946 // Triangular += Product
947 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
948 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::add_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
949 {
950  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
951  static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar,typename SrcXprType::Scalar> &)
952  {
953  dst._assignProduct(src, 1);
954  }
955 };
956 
957 // Triangular -= Product
958 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
959 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::sub_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
960 {
961  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
962  static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,typename SrcXprType::Scalar> &)
963  {
964  dst._assignProduct(src, -1);
965  }
966 };
967 
968 } // end namespace internal
969 
970 } // end namespace Eigen
971 
972 #endif // EIGEN_TRIANGULARMATRIX_H
const NestedExpression & nestedExpression() const
Definition: TriangularMatrix.h:233
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
Index innerStride() const
Definition: TriangularMatrix.h:364
const unsigned int LvalueBit
Definition: Constants.h:139
bool isLowerTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:668
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:271
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
Scalar determinant() const
Definition: TriangularMatrix.h:309
const unsigned int PacketAccessBit
Definition: Constants.h:89
Definition: EigenBase.h:28
const ConstTransposeReturnType transpose() const
Definition: TriangularMatrix.h:264
Scalar coeff(Index row, Index col) const
Definition: TriangularMatrix.h:406
const Product< TriangularViewType, OtherDerived > operator*(const MatrixBase< OtherDerived > &rhs) const
Definition: TriangularMatrix.h:453
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
const ConjugateReturnType conjugate() const
Definition: TriangularMatrix.h:242
Definition: Constants.h:208
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: XprHelper.h:35
Index rows() const
Definition: TriangularMatrix.h:226
void swap(TriangularBase< OtherDerived > &other)
Definition: TriangularMatrix.h:513
TriangularViewType & operator/=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:386
const AdjointReturnType adjoint() const
Definition: TriangularMatrix.h:248
Expression of a selfadjoint matrix from a triangular part of a dense matrix.
Definition: SelfAdjointView.h:49
Definition: Constants.h:218
void evalTo(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:595
Definition: Constants.h:214
Definition: Constants.h:206
Index cols() const
Definition: TriangularMatrix.h:229
TriangularViewType & setZero()
Definition: TriangularMatrix.h:397
Definition: Constants.h:210
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: Constants.h:212
Definition: Eigen_Colamd.h:50
Scalar & coeffRef(Index row, Index col)
Definition: TriangularMatrix.h:416
Definition: Constants.h:220
bool isUpperTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:642
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:186
const SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView() const
Definition: TriangularMatrix.h:299
SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView()
Definition: TriangularMatrix.h:291
TriangularViewType & operator-=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:376
Definition: Constants.h:491
const int Dynamic
Definition: Constants.h:21
Pseudo expression representing a solving operation.
Definition: Solve.h:62
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
TriangularViewType & operator+=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:369
void evalToLazy(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:926