Eigen  3.2.92
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>::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  // FIXME This, combined with const_cast_derived in transpose() leads to a const-correctness loophole
217  EIGEN_DEVICE_FUNC
218  explicit inline TriangularView(MatrixType& matrix) : m_matrix(matrix)
219  {}
220 
221  using Base::operator=;
222  TriangularView& operator=(const TriangularView &other)
223  { return Base::operator=(other); }
224 
226  EIGEN_DEVICE_FUNC
227  inline Index rows() const { return m_matrix.rows(); }
229  EIGEN_DEVICE_FUNC
230  inline Index cols() const { return m_matrix.cols(); }
231 
233  EIGEN_DEVICE_FUNC
234  const NestedExpression& nestedExpression() const { return m_matrix; }
235 
237  EIGEN_DEVICE_FUNC
238  NestedExpression& nestedExpression() { return *const_cast<NestedExpression*>(&m_matrix); }
239 
242  EIGEN_DEVICE_FUNC
243  inline const ConjugateReturnType conjugate() const
244  { return ConjugateReturnType(m_matrix.conjugate()); }
245 
248  EIGEN_DEVICE_FUNC
249  inline const AdjointReturnType adjoint() const
250  { return AdjointReturnType(m_matrix.adjoint()); }
251 
254  EIGEN_DEVICE_FUNC
255  inline TransposeReturnType transpose()
256  {
257  EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
258  typename MatrixType::TransposeReturnType tmp(m_matrix.const_cast_derived());
259  return TransposeReturnType(tmp);
260  }
261 
264  EIGEN_DEVICE_FUNC
265  inline const ConstTransposeReturnType transpose() const
266  {
267  return ConstTransposeReturnType(m_matrix.transpose());
268  }
269 
270  template<typename Other>
271  EIGEN_DEVICE_FUNC
272  inline const Solve<TriangularView, Other>
273  solve(const MatrixBase<Other>& other) const
274  { return Solve<TriangularView, Other>(*this, other.derived()); }
275 
276  // workaround MSVC ICE
277  #if EIGEN_COMP_MSVC
278  template<int Side, typename Other>
279  EIGEN_DEVICE_FUNC
280  inline const internal::triangular_solve_retval<Side,TriangularView, Other>
281  solve(const MatrixBase<Other>& other) const
282  { return Base::template solve<Side>(other); }
283  #else
284  using Base::solve;
285  #endif
286 
291  EIGEN_DEVICE_FUNC
293  {
294  EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
296  }
297 
299  EIGEN_DEVICE_FUNC
301  {
302  EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
304  }
305 
306 
309  EIGEN_DEVICE_FUNC
310  Scalar determinant() const
311  {
312  if (Mode & UnitDiag)
313  return 1;
314  else if (Mode & ZeroDiag)
315  return 0;
316  else
317  return m_matrix.diagonal().prod();
318  }
319 
320  protected:
321 
322  MatrixTypeNested m_matrix;
323 };
324 
334 template<typename _MatrixType, unsigned int _Mode> class TriangularViewImpl<_MatrixType,_Mode,Dense>
335  : public TriangularBase<TriangularView<_MatrixType, _Mode> >
336 {
337  public:
338 
341  typedef typename internal::traits<TriangularViewType>::Scalar Scalar;
342 
343  typedef _MatrixType MatrixType;
344  typedef typename MatrixType::PlainObject DenseMatrixType;
345  typedef DenseMatrixType PlainObject;
346 
347  public:
348  using Base::evalToLazy;
349  using Base::derived;
350 
351  typedef typename internal::traits<TriangularViewType>::StorageKind StorageKind;
352 
353  enum {
354  Mode = _Mode,
355  Flags = internal::traits<TriangularViewType>::Flags
356  };
357 
360  EIGEN_DEVICE_FUNC
361  inline Index outerStride() const { return derived().nestedExpression().outerStride(); }
364  EIGEN_DEVICE_FUNC
365  inline Index innerStride() const { return derived().nestedExpression().innerStride(); }
366 
368  template<typename Other>
369  EIGEN_DEVICE_FUNC
370  TriangularViewType& operator+=(const DenseBase<Other>& other) {
371  internal::call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op<Scalar>());
372  return derived();
373  }
375  template<typename Other>
376  EIGEN_DEVICE_FUNC
377  TriangularViewType& operator-=(const DenseBase<Other>& other) {
378  internal::call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op<Scalar>());
379  return derived();
380  }
381 
383  EIGEN_DEVICE_FUNC
384  TriangularViewType& operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() * other; }
386  EIGEN_DEVICE_FUNC
387  TriangularViewType& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() / other; }
388 
390  EIGEN_DEVICE_FUNC
391  void fill(const Scalar& value) { setConstant(value); }
393  EIGEN_DEVICE_FUNC
394  TriangularViewType& setConstant(const Scalar& value)
395  { return *this = MatrixType::Constant(derived().rows(), derived().cols(), value); }
397  EIGEN_DEVICE_FUNC
398  TriangularViewType& setZero() { return setConstant(Scalar(0)); }
400  EIGEN_DEVICE_FUNC
401  TriangularViewType& setOnes() { return setConstant(Scalar(1)); }
402 
406  EIGEN_DEVICE_FUNC
407  inline Scalar coeff(Index row, Index col) const
408  {
409  Base::check_coordinates_internal(row, col);
410  return derived().nestedExpression().coeff(row, col);
411  }
412 
416  EIGEN_DEVICE_FUNC
417  inline Scalar& coeffRef(Index row, Index col)
418  {
419  EIGEN_STATIC_ASSERT_LVALUE(TriangularViewType);
420  Base::check_coordinates_internal(row, col);
421  return derived().nestedExpression().const_cast_derived().coeffRef(row, col);
422  }
423 
425  template<typename OtherDerived>
426  EIGEN_DEVICE_FUNC
427  TriangularViewType& operator=(const TriangularBase<OtherDerived>& other);
428 
430  template<typename OtherDerived>
431  EIGEN_DEVICE_FUNC
432  TriangularViewType& operator=(const MatrixBase<OtherDerived>& other);
433 
434 #ifndef EIGEN_PARSED_BY_DOXYGEN
435  EIGEN_DEVICE_FUNC
436  TriangularViewType& operator=(const TriangularViewImpl& other)
437  { return *this = other.derived().nestedExpression(); }
438 
440  template<typename OtherDerived>
441  EIGEN_DEVICE_FUNC
442  void lazyAssign(const TriangularBase<OtherDerived>& other);
443 
445  template<typename OtherDerived>
446  EIGEN_DEVICE_FUNC
447  void lazyAssign(const MatrixBase<OtherDerived>& other);
448 #endif
449 
451  template<typename OtherDerived>
452  EIGEN_DEVICE_FUNC
455  {
456  return Product<TriangularViewType,OtherDerived>(derived(), rhs.derived());
457  }
458 
460  template<typename OtherDerived> friend
461  EIGEN_DEVICE_FUNC
463  operator*(const MatrixBase<OtherDerived>& lhs, const TriangularViewImpl& rhs)
464  {
465  return Product<OtherDerived,TriangularViewType>(lhs.derived(),rhs.derived());
466  }
467 
489  template<int Side, typename Other>
490  EIGEN_DEVICE_FUNC
491  inline const internal::triangular_solve_retval<Side,TriangularViewType, Other>
492  solve(const MatrixBase<Other>& other) const;
493 
501  template<int Side, typename OtherDerived>
502  EIGEN_DEVICE_FUNC
503  void solveInPlace(const MatrixBase<OtherDerived>& other) const;
504 
505  template<typename OtherDerived>
506  EIGEN_DEVICE_FUNC
507  void solveInPlace(const MatrixBase<OtherDerived>& other) const
508  { return solveInPlace<OnTheLeft>(other); }
509 
511  template<typename OtherDerived>
512  EIGEN_DEVICE_FUNC
513 #ifdef EIGEN_PARSED_BY_DOXYGEN
515 #else
516  void swap(TriangularBase<OtherDerived> const & other)
517 #endif
518  {
519  EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
520  call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
521  }
522 
525  template<typename OtherDerived>
526  EIGEN_DEVICE_FUNC
527  void swap(MatrixBase<OtherDerived> const & other)
528  {
529  EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
530  call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
531  }
532 
533  template<typename RhsType, typename DstType>
534  EIGEN_DEVICE_FUNC
535  EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const {
536  if(!(internal::is_same<RhsType,DstType>::value && internal::extract_data(dst) == internal::extract_data(rhs)))
537  dst = rhs;
538  this->solveInPlace(dst);
539  }
540 
541  template<typename ProductType>
542  EIGEN_DEVICE_FUNC
543  EIGEN_STRONG_INLINE TriangularViewType& _assignProduct(const ProductType& prod, const Scalar& alpha);
544 };
545 
546 /***************************************************************************
547 * Implementation of triangular evaluation/assignment
548 ***************************************************************************/
549 
550 // FIXME should we keep that possibility
551 template<typename MatrixType, unsigned int Mode>
552 template<typename OtherDerived>
554 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const MatrixBase<OtherDerived>& other)
555 {
556  internal::call_assignment_no_alias(derived(), other.derived(), internal::assign_op<Scalar>());
557  return derived();
558 }
559 
560 // FIXME should we keep that possibility
561 template<typename MatrixType, unsigned int Mode>
562 template<typename OtherDerived>
563 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const MatrixBase<OtherDerived>& other)
564 {
565  internal::call_assignment_no_alias(derived(), other.template triangularView<Mode>());
566 }
567 
568 
569 
570 template<typename MatrixType, unsigned int Mode>
571 template<typename OtherDerived>
573 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const TriangularBase<OtherDerived>& other)
574 {
575  eigen_assert(Mode == int(OtherDerived::Mode));
576  internal::call_assignment(derived(), other.derived());
577  return derived();
578 }
579 
580 template<typename MatrixType, unsigned int Mode>
581 template<typename OtherDerived>
582 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const TriangularBase<OtherDerived>& other)
583 {
584  eigen_assert(Mode == int(OtherDerived::Mode));
585  internal::call_assignment_no_alias(derived(), other.derived());
586 }
587 
588 /***************************************************************************
589 * Implementation of TriangularBase methods
590 ***************************************************************************/
591 
594 template<typename Derived>
595 template<typename DenseDerived>
597 {
598  if(internal::traits<Derived>::Flags & EvalBeforeAssigningBit)
599  {
600  typename internal::plain_matrix_type<Derived>::type other_evaluated(rows(), cols());
601  evalToLazy(other_evaluated);
602  other.derived().swap(other_evaluated);
603  }
604  else
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>
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  using std::abs;
653  RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
654  for(Index j = 0; j < cols(); ++j)
655  {
656  Index maxi = (std::min)(j, rows()-1);
657  for(Index i = 0; i <= maxi; ++i)
658  {
659  RealScalar absValue = abs(coeff(i,j));
660  if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
661  }
662  }
663  RealScalar threshold = maxAbsOnUpperPart * prec;
664  for(Index j = 0; j < cols(); ++j)
665  for(Index i = j+1; i < rows(); ++i)
666  if(abs(coeff(i, j)) > threshold) return false;
667  return true;
668 }
669 
675 template<typename Derived>
676 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const
677 {
678  using std::abs;
679  RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
680  for(Index j = 0; j < cols(); ++j)
681  for(Index i = j; i < rows(); ++i)
682  {
683  RealScalar absValue = abs(coeff(i,j));
684  if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
685  }
686  RealScalar threshold = maxAbsOnLowerPart * prec;
687  for(Index j = 1; j < cols(); ++j)
688  {
689  Index maxi = (std::min)(j, rows()-1);
690  for(Index i = 0; i < maxi; ++i)
691  if(abs(coeff(i, j)) > threshold) return false;
692  }
693  return true;
694 }
695 
696 
697 /***************************************************************************
698 ****************************************************************************
699 * Evaluators and Assignment of triangular expressions
700 ***************************************************************************
701 ***************************************************************************/
702 
703 namespace internal {
704 
705 
706 // TODO currently a triangular expression has the form TriangularView<.,.>
707 // in the future triangular-ness should be defined by the expression traits
708 // such that Transpose<TriangularView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
709 template<typename MatrixType, unsigned int Mode>
710 struct evaluator_traits<TriangularView<MatrixType,Mode> >
711 {
712  typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
713  typedef typename glue_shapes<typename evaluator_traits<MatrixType>::Shape, TriangularShape>::type Shape;
714 };
715 
716 template<typename MatrixType, unsigned int Mode>
717 struct unary_evaluator<TriangularView<MatrixType,Mode>, IndexBased>
718  : evaluator<typename internal::remove_all<MatrixType>::type>
719 {
720  typedef TriangularView<MatrixType,Mode> XprType;
721  typedef evaluator<typename internal::remove_all<MatrixType>::type> Base;
722  unary_evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {}
723 };
724 
725 // Additional assignment kinds:
726 struct Triangular2Triangular {};
727 struct Triangular2Dense {};
728 struct Dense2Triangular {};
729 
730 
731 template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop;
732 
733 
739 template<int UpLo, int Mode, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version = Specialized>
740 class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version>
741 {
742 protected:
743  typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base;
744  typedef typename Base::DstXprType DstXprType;
745  typedef typename Base::SrcXprType SrcXprType;
746  using Base::m_dst;
747  using Base::m_src;
748  using Base::m_functor;
749 public:
750 
751  typedef typename Base::DstEvaluatorType DstEvaluatorType;
752  typedef typename Base::SrcEvaluatorType SrcEvaluatorType;
753  typedef typename Base::Scalar Scalar;
754  typedef typename Base::AssignmentTraits AssignmentTraits;
755 
756 
757  EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr)
758  : Base(dst, src, func, dstExpr)
759  {}
760 
761 #ifdef EIGEN_INTERNAL_DEBUGGING
762  EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col)
763  {
764  eigen_internal_assert(row!=col);
765  Base::assignCoeff(row,col);
766  }
767 #else
768  using Base::assignCoeff;
769 #endif
770 
771  EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id)
772  {
773  if(Mode==UnitDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(1));
774  else if(Mode==ZeroDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(0));
775  else if(Mode==0) Base::assignCoeff(id,id);
776  }
777 
778  EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index row, Index col)
779  {
780  eigen_internal_assert(row!=col);
781  if(SetOpposite)
782  m_functor.assignCoeff(m_dst.coeffRef(row,col), Scalar(0));
783  }
784 };
785 
786 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType, typename Functor>
787 EIGEN_DEVICE_FUNC void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src, const Functor &func)
788 {
789  eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
790 
791  typedef evaluator<DstXprType> DstEvaluatorType;
792  typedef evaluator<SrcXprType> SrcEvaluatorType;
793 
794  DstEvaluatorType dstEvaluator(dst);
795  SrcEvaluatorType srcEvaluator(src);
796 
797  typedef triangular_dense_assignment_kernel< Mode&(Lower|Upper),Mode&(UnitDiag|ZeroDiag|SelfAdjoint),SetOpposite,
798  DstEvaluatorType,SrcEvaluatorType,Functor> Kernel;
799  Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived());
800 
801  enum {
802  unroll = DstXprType::SizeAtCompileTime != Dynamic
803  && SrcEvaluatorType::CoeffReadCost < HugeCost
804  && DstXprType::SizeAtCompileTime * SrcEvaluatorType::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT
805  };
806 
807  triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, SetOpposite>::run(kernel);
808 }
809 
810 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType>
811 EIGEN_DEVICE_FUNC void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src)
812 {
813  call_triangular_assignment_loop<Mode,SetOpposite>(dst, src, internal::assign_op<typename DstXprType::Scalar>());
814 }
815 
816 template<> struct AssignmentKind<TriangularShape,TriangularShape> { typedef Triangular2Triangular Kind; };
817 template<> struct AssignmentKind<DenseShape,TriangularShape> { typedef Triangular2Dense Kind; };
818 template<> struct AssignmentKind<TriangularShape,DenseShape> { typedef Dense2Triangular Kind; };
819 
820 
821 template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar>
822 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Triangular, Scalar>
823 {
824  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
825  {
826  eigen_assert(int(DstXprType::Mode) == int(SrcXprType::Mode));
827 
828  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
829  }
830 };
831 
832 template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar>
833 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Dense, Scalar>
834 {
835  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
836  {
837  call_triangular_assignment_loop<SrcXprType::Mode, (SrcXprType::Mode&SelfAdjoint)==0>(dst, src, func);
838  }
839 };
840 
841 template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar>
842 struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular, Scalar>
843 {
844  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
845  {
846  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
847  }
848 };
849 
850 
851 template<typename Kernel, unsigned int Mode, int UnrollCount, bool SetOpposite>
852 struct triangular_assignment_loop
853 {
854  // FIXME: this is not very clean, perhaps this information should be provided by the kernel?
855  typedef typename Kernel::DstEvaluatorType DstEvaluatorType;
856  typedef typename DstEvaluatorType::XprType DstXprType;
857 
858  enum {
859  col = (UnrollCount-1) / DstXprType::RowsAtCompileTime,
860  row = (UnrollCount-1) % DstXprType::RowsAtCompileTime
861  };
862 
863  typedef typename Kernel::Scalar Scalar;
864 
865  EIGEN_DEVICE_FUNC
866  static inline void run(Kernel &kernel)
867  {
868  triangular_assignment_loop<Kernel, Mode, UnrollCount-1, SetOpposite>::run(kernel);
869 
870  if(row==col)
871  kernel.assignDiagonalCoeff(row);
872  else if( ((Mode&Lower) && row>col) || ((Mode&Upper) && row<col) )
873  kernel.assignCoeff(row,col);
874  else if(SetOpposite)
875  kernel.assignOppositeCoeff(row,col);
876  }
877 };
878 
879 // prevent buggy user code from causing an infinite recursion
880 template<typename Kernel, unsigned int Mode, bool SetOpposite>
881 struct triangular_assignment_loop<Kernel, Mode, 0, SetOpposite>
882 {
883  EIGEN_DEVICE_FUNC
884  static inline void run(Kernel &) {}
885 };
886 
887 
888 
889 // TODO: experiment with a recursive assignment procedure splitting the current
890 // triangular part into one rectangular and two triangular parts.
891 
892 
893 template<typename Kernel, unsigned int Mode, bool SetOpposite>
894 struct triangular_assignment_loop<Kernel, Mode, Dynamic, SetOpposite>
895 {
896  typedef typename Kernel::Scalar Scalar;
897  EIGEN_DEVICE_FUNC
898  static inline void run(Kernel &kernel)
899  {
900  for(Index j = 0; j < kernel.cols(); ++j)
901  {
902  Index maxi = (std::min)(j, kernel.rows());
903  Index i = 0;
904  if (((Mode&Lower) && SetOpposite) || (Mode&Upper))
905  {
906  for(; i < maxi; ++i)
907  if(Mode&Upper) kernel.assignCoeff(i, j);
908  else kernel.assignOppositeCoeff(i, j);
909  }
910  else
911  i = maxi;
912 
913  if(i<kernel.rows()) // then i==j
914  kernel.assignDiagonalCoeff(i++);
915 
916  if (((Mode&Upper) && SetOpposite) || (Mode&Lower))
917  {
918  for(; i < kernel.rows(); ++i)
919  if(Mode&Lower) kernel.assignCoeff(i, j);
920  else kernel.assignOppositeCoeff(i, j);
921  }
922  }
923  }
924 };
925 
926 } // end namespace internal
927 
930 template<typename Derived>
931 template<typename DenseDerived>
933 {
934  other.derived().resize(this->rows(), this->cols());
935  internal::call_triangular_assignment_loop<Derived::Mode,(Derived::Mode&SelfAdjoint)==0 /* SetOpposite */>(other.derived(), derived().nestedExpression());
936 }
937 
938 namespace internal {
939 
940 // Triangular = Product
941 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
942 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::assign_op<Scalar>, Dense2Triangular, Scalar>
943 {
944  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
945  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
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>, Dense2Triangular, Scalar>
955 {
956  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
957  static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<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>, Dense2Triangular, Scalar>
966 {
967  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
968  static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<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 NestedExpression & nestedExpression() const
Definition: TriangularMatrix.h:234
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:107
Index outerStride() const
Definition: TriangularMatrix.h:361
NestedExpression & nestedExpression()
Definition: TriangularMatrix.h:238
TriangularViewType & setConstant(const Scalar &value)
Definition: TriangularMatrix.h:394
TransposeReturnType transpose()
Definition: TriangularMatrix.h:255
Base class for triangular part in a matrix.
Definition: TriangularMatrix.h:27
const unsigned int DirectAccessBit
Definition: Constants.h:149
Index innerStride() const
Definition: TriangularMatrix.h:365
const unsigned int LvalueBit
Definition: Constants.h:138
bool isLowerTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:676
Definition: Constants.h:210
friend const Product< OtherDerived, TriangularViewType > operator*(const MatrixBase< OtherDerived > &lhs, const TriangularViewImpl &rhs)
Definition: TriangularMatrix.h:463
Definition: LDLT.h:16
TriangularViewType & setOnes()
Definition: TriangularMatrix.h:401
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:384
Base class for all dense matrices, vectors, and arrays.
Definition: DenseBase.h:41
Scalar determinant() const
Definition: TriangularMatrix.h:310
const unsigned int PacketAccessBit
Definition: Constants.h:88
void resize(Index newSize)
Definition: DenseBase.h:245
Definition: EigenBase.h:28
const ConstTransposeReturnType transpose() const
Definition: TriangularMatrix.h:265
Scalar coeff(Index row, Index col) const
Definition: TriangularMatrix.h:407
const Product< TriangularViewType, OtherDerived > operator*(const MatrixBase< OtherDerived > &rhs) const
Definition: TriangularMatrix.h:454
Definition: Constants.h:214
Definition: Constants.h:204
void copyCoeff(Index row, Index col, Other &other)
Definition: TriangularMatrix.h:84
void swap(MatrixBase< OtherDerived > const &other)
Definition: TriangularMatrix.h:527
Definition: Constants.h:212
Definition: Constants.h:218
const ConjugateReturnType conjugate() const
Definition: TriangularMatrix.h:243
Index rows() const
Definition: TriangularMatrix.h:227
void swap(TriangularBase< OtherDerived > &other)
Definition: TriangularMatrix.h:514
TriangularViewType & operator/=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:387
const AdjointReturnType adjoint() const
Definition: TriangularMatrix.h:249
Expression of a selfadjoint matrix from a triangular part of a dense matrix.
Definition: SelfAdjointView.h:49
const unsigned int EvalBeforeAssigningBit
Definition: Constants.h:70
void evalTo(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:596
void swap(const DenseBase< OtherDerived > &other)
Definition: DenseBase.h:418
Index cols() const
Definition: TriangularMatrix.h:230
TriangularViewType & setZero()
Definition: TriangularMatrix.h:398
Definition: Eigen_Colamd.h:54
Scalar & coeffRef(Index row, Index col)
Definition: TriangularMatrix.h:417
bool isUpperTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:650
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:186
const SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView() const
Definition: TriangularMatrix.h:300
SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView()
Definition: TriangularMatrix.h:292
Definition: Constants.h:208
TriangularViewType & operator-=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:377
Definition: Constants.h:490
Definition: Constants.h:216
Pseudo expression representing a solving operation.
Definition: Solve.h:62
Definition: Constants.h:220
Definition: Constants.h:206
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const unsigned int LinearAccessBit
Definition: Constants.h:124
void fill(const Scalar &value)
Definition: TriangularMatrix.h:391
TriangularViewType & operator+=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:370
void evalToLazy(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:932