Eigen  3.2.91
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 
225  EIGEN_DEVICE_FUNC
226  inline Index rows() const { return m_matrix.rows(); }
227  EIGEN_DEVICE_FUNC
228  inline Index cols() const { return m_matrix.cols(); }
229 
230  EIGEN_DEVICE_FUNC
231  const NestedExpression& nestedExpression() const { return m_matrix; }
232  EIGEN_DEVICE_FUNC
233  NestedExpression& nestedExpression() { return *const_cast<NestedExpression*>(&m_matrix); }
234 
237  EIGEN_DEVICE_FUNC
238  inline const ConjugateReturnType conjugate() const
239  { return ConjugateReturnType(m_matrix.conjugate()); }
240 
243  EIGEN_DEVICE_FUNC
244  inline const AdjointReturnType adjoint() const
245  { return AdjointReturnType(m_matrix.adjoint()); }
246 
249  EIGEN_DEVICE_FUNC
250  inline TransposeReturnType transpose()
251  {
252  EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
253  typename MatrixType::TransposeReturnType tmp(m_matrix.const_cast_derived());
254  return TransposeReturnType(tmp);
255  }
256 
259  EIGEN_DEVICE_FUNC
260  inline const ConstTransposeReturnType transpose() const
261  {
262  return ConstTransposeReturnType(m_matrix.transpose());
263  }
264 
265  template<typename Other>
266  EIGEN_DEVICE_FUNC
267  inline const Solve<TriangularView, Other>
268  solve(const MatrixBase<Other>& other) const
269  { return Solve<TriangularView, Other>(*this, other.derived()); }
270 
271  // workaround MSVC ICE
272  #if EIGEN_COMP_MSVC
273  template<int Side, typename Other>
274  EIGEN_DEVICE_FUNC
275  inline const internal::triangular_solve_retval<Side,TriangularView, Other>
276  solve(const MatrixBase<Other>& other) const
277  { return Base::template solve<Side>(other); }
278  #else
279  using Base::solve;
280  #endif
281 
282  EIGEN_DEVICE_FUNC
283  const SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() const
284  {
285  EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
286  return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
287  }
288  EIGEN_DEVICE_FUNC
289  SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView()
290  {
291  EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
292  return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
293  }
294 
295  EIGEN_DEVICE_FUNC
296  Scalar determinant() const
297  {
298  if (Mode & UnitDiag)
299  return 1;
300  else if (Mode & ZeroDiag)
301  return 0;
302  else
303  return m_matrix.diagonal().prod();
304  }
305 
306  protected:
307 
308  MatrixTypeNested m_matrix;
309 };
310 
320 template<typename _MatrixType, unsigned int _Mode> class TriangularViewImpl<_MatrixType,_Mode,Dense>
321  : public TriangularBase<TriangularView<_MatrixType, _Mode> >
322 {
323  public:
324 
327  typedef typename internal::traits<TriangularViewType>::Scalar Scalar;
328 
329  typedef _MatrixType MatrixType;
330  typedef typename MatrixType::PlainObject DenseMatrixType;
331  typedef DenseMatrixType PlainObject;
332 
333  public:
334  using Base::evalToLazy;
335  using Base::derived;
336 
337  typedef typename internal::traits<TriangularViewType>::StorageKind StorageKind;
338 
339  enum {
340  Mode = _Mode,
341  Flags = internal::traits<TriangularViewType>::Flags
342  };
343 
344  EIGEN_DEVICE_FUNC
345  inline Index outerStride() const { return derived().nestedExpression().outerStride(); }
346  EIGEN_DEVICE_FUNC
347  inline Index innerStride() const { return derived().nestedExpression().innerStride(); }
348 
350  template<typename Other>
351  EIGEN_DEVICE_FUNC
352  TriangularViewType& operator+=(const DenseBase<Other>& other) {
353  internal::call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op<Scalar>());
354  return derived();
355  }
357  template<typename Other>
358  EIGEN_DEVICE_FUNC
359  TriangularViewType& operator-=(const DenseBase<Other>& other) {
360  internal::call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op<Scalar>());
361  return derived();
362  }
363 
365  EIGEN_DEVICE_FUNC
366  TriangularViewType& operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() * other; }
368  EIGEN_DEVICE_FUNC
369  TriangularViewType& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() / other; }
370 
372  EIGEN_DEVICE_FUNC
373  void fill(const Scalar& value) { setConstant(value); }
375  EIGEN_DEVICE_FUNC
376  TriangularViewType& setConstant(const Scalar& value)
377  { return *this = MatrixType::Constant(derived().rows(), derived().cols(), value); }
379  EIGEN_DEVICE_FUNC
380  TriangularViewType& setZero() { return setConstant(Scalar(0)); }
382  EIGEN_DEVICE_FUNC
383  TriangularViewType& setOnes() { return setConstant(Scalar(1)); }
384 
388  EIGEN_DEVICE_FUNC
389  inline Scalar coeff(Index row, Index col) const
390  {
391  Base::check_coordinates_internal(row, col);
392  return derived().nestedExpression().coeff(row, col);
393  }
394 
398  EIGEN_DEVICE_FUNC
399  inline Scalar& coeffRef(Index row, Index col)
400  {
401  EIGEN_STATIC_ASSERT_LVALUE(TriangularViewType);
402  Base::check_coordinates_internal(row, col);
403  return derived().nestedExpression().const_cast_derived().coeffRef(row, col);
404  }
405 
407  template<typename OtherDerived>
408  EIGEN_DEVICE_FUNC
409  TriangularViewType& operator=(const TriangularBase<OtherDerived>& other);
410 
411  template<typename OtherDerived>
412  EIGEN_DEVICE_FUNC
413  TriangularViewType& operator=(const MatrixBase<OtherDerived>& other);
414 
415  EIGEN_DEVICE_FUNC
416  TriangularViewType& operator=(const TriangularViewImpl& other)
417  { return *this = other.derived().nestedExpression(); }
418 
419  template<typename OtherDerived>
420  EIGEN_DEVICE_FUNC
421  void lazyAssign(const TriangularBase<OtherDerived>& other);
422 
423  template<typename OtherDerived>
424  EIGEN_DEVICE_FUNC
425  void lazyAssign(const MatrixBase<OtherDerived>& other);
426 
428  template<typename OtherDerived>
429  EIGEN_DEVICE_FUNC
430  const Product<TriangularViewType,OtherDerived>
432  {
433  return Product<TriangularViewType,OtherDerived>(derived(), rhs.derived());
434  }
435 
437  template<typename OtherDerived> friend
438  EIGEN_DEVICE_FUNC
440  operator*(const MatrixBase<OtherDerived>& lhs, const TriangularViewImpl& rhs)
441  {
442  return Product<OtherDerived,TriangularViewType>(lhs.derived(),rhs.derived());
443  }
444 
445  template<int Side, typename Other>
446  EIGEN_DEVICE_FUNC
447  inline const internal::triangular_solve_retval<Side,TriangularViewType, Other>
448  solve(const MatrixBase<Other>& other) const;
449 
450  template<int Side, typename OtherDerived>
451  EIGEN_DEVICE_FUNC
452  void solveInPlace(const MatrixBase<OtherDerived>& other) const;
453 
454  template<typename OtherDerived>
455  EIGEN_DEVICE_FUNC
456  void solveInPlace(const MatrixBase<OtherDerived>& other) const
457  { return solveInPlace<OnTheLeft>(other); }
458 
459  template<typename OtherDerived>
460  EIGEN_DEVICE_FUNC
461  void swap(TriangularBase<OtherDerived> const & other)
462  {
463  call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
464  }
465 
466  // TODO: this overload is ambiguous and it should be deprecated (Gael)
467  template<typename OtherDerived>
468  EIGEN_DEVICE_FUNC
469  void swap(MatrixBase<OtherDerived> const & other)
470  {
471  call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
472  }
473 
474  template<typename RhsType, typename DstType>
475  EIGEN_DEVICE_FUNC
476  EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const {
477  if(!(internal::is_same<RhsType,DstType>::value && internal::extract_data(dst) == internal::extract_data(rhs)))
478  dst = rhs;
479  this->solveInPlace(dst);
480  }
481 
482  template<typename ProductType>
483  EIGEN_DEVICE_FUNC
484  EIGEN_STRONG_INLINE TriangularViewType& _assignProduct(const ProductType& prod, const Scalar& alpha);
485 };
486 
487 /***************************************************************************
488 * Implementation of triangular evaluation/assignment
489 ***************************************************************************/
490 
491 // FIXME should we keep that possibility
492 template<typename MatrixType, unsigned int Mode>
493 template<typename OtherDerived>
494 inline TriangularView<MatrixType, Mode>&
495 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const MatrixBase<OtherDerived>& other)
496 {
497  internal::call_assignment_no_alias(derived(), other.derived(), internal::assign_op<Scalar>());
498  return derived();
499 }
500 
501 // FIXME should we keep that possibility
502 template<typename MatrixType, unsigned int Mode>
503 template<typename OtherDerived>
504 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const MatrixBase<OtherDerived>& other)
505 {
506  internal::call_assignment(derived().noalias(), other.template triangularView<Mode>());
507 }
508 
509 
510 
511 template<typename MatrixType, unsigned int Mode>
512 template<typename OtherDerived>
513 inline TriangularView<MatrixType, Mode>&
514 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const TriangularBase<OtherDerived>& other)
515 {
516  eigen_assert(Mode == int(OtherDerived::Mode));
517  internal::call_assignment(derived(), other.derived());
518  return derived();
519 }
520 
521 template<typename MatrixType, unsigned int Mode>
522 template<typename OtherDerived>
523 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const TriangularBase<OtherDerived>& other)
524 {
525  eigen_assert(Mode == int(OtherDerived::Mode));
526  internal::call_assignment(derived().noalias(), other.derived());
527 }
528 
529 /***************************************************************************
530 * Implementation of TriangularBase methods
531 ***************************************************************************/
532 
535 template<typename Derived>
536 template<typename DenseDerived>
538 {
539  if(internal::traits<Derived>::Flags & EvalBeforeAssigningBit)
540  {
541  typename internal::plain_matrix_type<Derived>::type other_evaluated(rows(), cols());
542  evalToLazy(other_evaluated);
543  other.derived().swap(other_evaluated);
544  }
545  else
546  evalToLazy(other.derived());
547 }
548 
549 /***************************************************************************
550 * Implementation of TriangularView methods
551 ***************************************************************************/
552 
553 /***************************************************************************
554 * Implementation of MatrixBase methods
555 ***************************************************************************/
556 
568 template<typename Derived>
569 template<unsigned int Mode>
570 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
572 {
573  return typename TriangularViewReturnType<Mode>::Type(derived());
574 }
575 
577 template<typename Derived>
578 template<unsigned int Mode>
579 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
581 {
582  return typename ConstTriangularViewReturnType<Mode>::Type(derived());
583 }
584 
590 template<typename Derived>
591 bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const
592 {
593  using std::abs;
594  RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
595  for(Index j = 0; j < cols(); ++j)
596  {
597  Index maxi = (std::min)(j, rows()-1);
598  for(Index i = 0; i <= maxi; ++i)
599  {
600  RealScalar absValue = abs(coeff(i,j));
601  if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
602  }
603  }
604  RealScalar threshold = maxAbsOnUpperPart * prec;
605  for(Index j = 0; j < cols(); ++j)
606  for(Index i = j+1; i < rows(); ++i)
607  if(abs(coeff(i, j)) > threshold) return false;
608  return true;
609 }
610 
616 template<typename Derived>
617 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const
618 {
619  using std::abs;
620  RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
621  for(Index j = 0; j < cols(); ++j)
622  for(Index i = j; i < rows(); ++i)
623  {
624  RealScalar absValue = abs(coeff(i,j));
625  if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
626  }
627  RealScalar threshold = maxAbsOnLowerPart * prec;
628  for(Index j = 1; j < cols(); ++j)
629  {
630  Index maxi = (std::min)(j, rows()-1);
631  for(Index i = 0; i < maxi; ++i)
632  if(abs(coeff(i, j)) > threshold) return false;
633  }
634  return true;
635 }
636 
637 
638 /***************************************************************************
639 ****************************************************************************
640 * Evaluators and Assignment of triangular expressions
641 ***************************************************************************
642 ***************************************************************************/
643 
644 namespace internal {
645 
646 
647 // TODO currently a triangular expression has the form TriangularView<.,.>
648 // in the future triangular-ness should be defined by the expression traits
649 // such that Transpose<TriangularView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
650 template<typename MatrixType, unsigned int Mode>
651 struct evaluator_traits<TriangularView<MatrixType,Mode> >
652 {
653  typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
654  typedef typename glue_shapes<typename evaluator_traits<MatrixType>::Shape, TriangularShape>::type Shape;
655 
656  // 1 if assignment A = B assumes aliasing when B is of type T and thus B needs to be evaluated into a
657  // temporary; 0 if not.
658  static const int AssumeAliasing = 0;
659 };
660 
661 template<typename MatrixType, unsigned int Mode>
662 struct unary_evaluator<TriangularView<MatrixType,Mode>, IndexBased>
663  : evaluator<typename internal::remove_all<MatrixType>::type>
664 {
665  typedef TriangularView<MatrixType,Mode> XprType;
666  typedef evaluator<typename internal::remove_all<MatrixType>::type> Base;
667  unary_evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {}
668 };
669 
670 // Additional assignment kinds:
671 struct Triangular2Triangular {};
672 struct Triangular2Dense {};
673 struct Dense2Triangular {};
674 
675 
676 template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop;
677 
678 
684 template<int UpLo, int Mode, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version = Specialized>
685 class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version>
686 {
687 protected:
688  typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base;
689  typedef typename Base::DstXprType DstXprType;
690  typedef typename Base::SrcXprType SrcXprType;
691  using Base::m_dst;
692  using Base::m_src;
693  using Base::m_functor;
694 public:
695 
696  typedef typename Base::DstEvaluatorType DstEvaluatorType;
697  typedef typename Base::SrcEvaluatorType SrcEvaluatorType;
698  typedef typename Base::Scalar Scalar;
699  typedef typename Base::AssignmentTraits AssignmentTraits;
700 
701 
702  EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr)
703  : Base(dst, src, func, dstExpr)
704  {}
705 
706 #ifdef EIGEN_INTERNAL_DEBUGGING
707  EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col)
708  {
709  eigen_internal_assert(row!=col);
710  Base::assignCoeff(row,col);
711  }
712 #else
713  using Base::assignCoeff;
714 #endif
715 
716  EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id)
717  {
718  if(Mode==UnitDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(1));
719  else if(Mode==ZeroDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(0));
720  else if(Mode==0) Base::assignCoeff(id,id);
721  }
722 
723  EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index row, Index col)
724  {
725  eigen_internal_assert(row!=col);
726  if(SetOpposite)
727  m_functor.assignCoeff(m_dst.coeffRef(row,col), Scalar(0));
728  }
729 };
730 
731 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType, typename Functor>
732 EIGEN_DEVICE_FUNC void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src, const Functor &func)
733 {
734  eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
735 
736  typedef evaluator<DstXprType> DstEvaluatorType;
737  typedef evaluator<SrcXprType> SrcEvaluatorType;
738 
739  DstEvaluatorType dstEvaluator(dst);
740  SrcEvaluatorType srcEvaluator(src);
741 
742  typedef triangular_dense_assignment_kernel< Mode&(Lower|Upper),Mode&(UnitDiag|ZeroDiag|SelfAdjoint),SetOpposite,
743  DstEvaluatorType,SrcEvaluatorType,Functor> Kernel;
744  Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived());
745 
746  enum {
747  unroll = DstXprType::SizeAtCompileTime != Dynamic
748  && SrcEvaluatorType::CoeffReadCost != Dynamic
749  && DstXprType::SizeAtCompileTime * SrcEvaluatorType::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT
750  };
751 
752  triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, SetOpposite>::run(kernel);
753 }
754 
755 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType>
756 EIGEN_DEVICE_FUNC void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src)
757 {
758  call_triangular_assignment_loop<Mode,SetOpposite>(dst, src, internal::assign_op<typename DstXprType::Scalar>());
759 }
760 
761 template<> struct AssignmentKind<TriangularShape,TriangularShape> { typedef Triangular2Triangular Kind; };
762 template<> struct AssignmentKind<DenseShape,TriangularShape> { typedef Triangular2Dense Kind; };
763 template<> struct AssignmentKind<TriangularShape,DenseShape> { typedef Dense2Triangular Kind; };
764 
765 
766 template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar>
767 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Triangular, Scalar>
768 {
769  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
770  {
771  eigen_assert(int(DstXprType::Mode) == int(SrcXprType::Mode));
772 
773  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
774  }
775 };
776 
777 template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar>
778 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Dense, Scalar>
779 {
780  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
781  {
782  call_triangular_assignment_loop<SrcXprType::Mode, (SrcXprType::Mode&SelfAdjoint)==0>(dst, src, func);
783  }
784 };
785 
786 template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar>
787 struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular, Scalar>
788 {
789  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
790  {
791  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
792  }
793 };
794 
795 
796 template<typename Kernel, unsigned int Mode, int UnrollCount, bool SetOpposite>
797 struct triangular_assignment_loop
798 {
799  // FIXME: this is not very clean, perhaps this information should be provided by the kernel?
800  typedef typename Kernel::DstEvaluatorType DstEvaluatorType;
801  typedef typename DstEvaluatorType::XprType DstXprType;
802 
803  enum {
804  col = (UnrollCount-1) / DstXprType::RowsAtCompileTime,
805  row = (UnrollCount-1) % DstXprType::RowsAtCompileTime
806  };
807 
808  typedef typename Kernel::Scalar Scalar;
809 
810  EIGEN_DEVICE_FUNC
811  static inline void run(Kernel &kernel)
812  {
813  triangular_assignment_loop<Kernel, Mode, UnrollCount-1, SetOpposite>::run(kernel);
814 
815  if(row==col)
816  kernel.assignDiagonalCoeff(row);
817  else if( ((Mode&Lower) && row>col) || ((Mode&Upper) && row<col) )
818  kernel.assignCoeff(row,col);
819  else if(SetOpposite)
820  kernel.assignOppositeCoeff(row,col);
821  }
822 };
823 
824 // prevent buggy user code from causing an infinite recursion
825 template<typename Kernel, unsigned int Mode, bool SetOpposite>
826 struct triangular_assignment_loop<Kernel, Mode, 0, SetOpposite>
827 {
828  EIGEN_DEVICE_FUNC
829  static inline void run(Kernel &) {}
830 };
831 
832 
833 
834 // TODO: experiment with a recursive assignment procedure splitting the current
835 // triangular part into one rectangular and two triangular parts.
836 
837 
838 template<typename Kernel, unsigned int Mode, bool SetOpposite>
839 struct triangular_assignment_loop<Kernel, Mode, Dynamic, SetOpposite>
840 {
841  typedef typename Kernel::Scalar Scalar;
842  EIGEN_DEVICE_FUNC
843  static inline void run(Kernel &kernel)
844  {
845  for(Index j = 0; j < kernel.cols(); ++j)
846  {
847  Index maxi = (std::min)(j, kernel.rows());
848  Index i = 0;
849  if (((Mode&Lower) && SetOpposite) || (Mode&Upper))
850  {
851  for(; i < maxi; ++i)
852  if(Mode&Upper) kernel.assignCoeff(i, j);
853  else kernel.assignOppositeCoeff(i, j);
854  }
855  else
856  i = maxi;
857 
858  if(i<kernel.rows()) // then i==j
859  kernel.assignDiagonalCoeff(i++);
860 
861  if (((Mode&Upper) && SetOpposite) || (Mode&Lower))
862  {
863  for(; i < kernel.rows(); ++i)
864  if(Mode&Lower) kernel.assignCoeff(i, j);
865  else kernel.assignOppositeCoeff(i, j);
866  }
867  }
868  }
869 };
870 
871 } // end namespace internal
872 
875 template<typename Derived>
876 template<typename DenseDerived>
878 {
879  other.derived().resize(this->rows(), this->cols());
880  internal::call_triangular_assignment_loop<Derived::Mode,(Derived::Mode&SelfAdjoint)==0 /* SetOpposite */>(other.derived(), derived().nestedExpression());
881 }
882 
883 namespace internal {
884 
885 // Triangular = Product
886 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
887 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::assign_op<Scalar>, Dense2Triangular, Scalar>
888 {
889  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
890  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
891  {
892  dst.setZero();
893  dst._assignProduct(src, 1);
894  }
895 };
896 
897 // Triangular += Product
898 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
899 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::add_assign_op<Scalar>, Dense2Triangular, Scalar>
900 {
901  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
902  static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar> &)
903  {
904  dst._assignProduct(src, 1);
905  }
906 };
907 
908 // Triangular -= Product
909 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
910 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::sub_assign_op<Scalar>, Dense2Triangular, Scalar>
911 {
912  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
913  static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar> &)
914  {
915  dst._assignProduct(src, -1);
916  }
917 };
918 
919 } // end namespace internal
920 
921 } // end namespace Eigen
922 
923 #endif // EIGEN_TRIANGULARMATRIX_H
Definition: Constants.h:204
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:107
TriangularViewType & setConstant(const Scalar &value)
Definition: TriangularMatrix.h:376
TransposeReturnType transpose()
Definition: TriangularMatrix.h:250
Base class for triangular part in a matrix.
Definition: TriangularMatrix.h:27
const unsigned int DirectAccessBit
Definition: Constants.h:141
Definition: Constants.h:196
Definition: Constants.h:202
const unsigned int LvalueBit
Definition: Constants.h:130
bool isLowerTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:617
friend const Product< OtherDerived, TriangularViewType > operator*(const MatrixBase< OtherDerived > &lhs, const TriangularViewImpl &rhs)
Definition: TriangularMatrix.h:440
Definition: LDLT.h:16
TriangularViewType & setOnes()
Definition: TriangularMatrix.h:383
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:366
Base class for all dense matrices, vectors, and arrays.
Definition: DenseBase.h:41
const unsigned int PacketAccessBit
Definition: Constants.h:80
void resize(Index newSize)
Definition: DenseBase.h:252
Definition: EigenBase.h:28
const ConstTransposeReturnType transpose() const
Definition: TriangularMatrix.h:260
Scalar coeff(Index row, Index col) const
Definition: TriangularMatrix.h:389
const Product< TriangularViewType, OtherDerived > operator*(const MatrixBase< OtherDerived > &rhs) const
Definition: TriangularMatrix.h:431
Definition: Constants.h:198
void copyCoeff(Index row, Index col, Other &other)
Definition: TriangularMatrix.h:84
TriangularView< const MatrixConjugateReturnType, Mode > ConjugateReturnType
Definition: TriangularMatrix.h:236
TriangularViewType & operator/=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:369
const AdjointReturnType adjoint() const
Definition: TriangularMatrix.h:244
Definition: Constants.h:210
const unsigned int EvalBeforeAssigningBit
Definition: Constants.h:62
void evalTo(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:537
void swap(const DenseBase< OtherDerived > &other)
Definition: DenseBase.h:425
Definition: Constants.h:200
TriangularViewType & setZero()
Definition: TriangularMatrix.h:380
Definition: Eigen_Colamd.h:54
Definition: Constants.h:208
Scalar & coeffRef(Index row, Index col)
Definition: TriangularMatrix.h:399
bool isUpperTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:591
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:186
Definition: Constants.h:206
TriangularViewType & operator-=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:359
Definition: Constants.h:482
Pseudo expression representing a solving operation.
Definition: Solve.h:63
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Definition: TriangularMatrix.h:38
const unsigned int LinearAccessBit
Definition: Constants.h:116
void fill(const Scalar &value)
Definition: TriangularMatrix.h:373
Definition: Constants.h:212
TriangularViewType & operator+=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:352
void evalToLazy(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:877