Eigen  3.3.7
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 
239  typedef TriangularView<const MatrixConjugateReturnType,Mode> ConjugateReturnType;
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
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
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
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
452  const Product<TriangularViewType,OtherDerived>
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 
490  template<int Side, typename Other>
491  EIGEN_DEVICE_FUNC
492  inline const internal::triangular_solve_retval<Side,TriangularViewType, Other>
493  solve(const MatrixBase<Other>& other) const;
494 
504  template<int Side, typename OtherDerived>
505  EIGEN_DEVICE_FUNC
506  void solveInPlace(const MatrixBase<OtherDerived>& other) const;
507 
508  template<typename OtherDerived>
509  EIGEN_DEVICE_FUNC
510  void solveInPlace(const MatrixBase<OtherDerived>& other) const
511  { return solveInPlace<OnTheLeft>(other); }
512 
514  template<typename OtherDerived>
515  EIGEN_DEVICE_FUNC
516 #ifdef EIGEN_PARSED_BY_DOXYGEN
518 #else
519  void swap(TriangularBase<OtherDerived> const & other)
520 #endif
521  {
522  EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
523  call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
524  }
525 
528  template<typename OtherDerived>
529  EIGEN_DEVICE_FUNC
530  void swap(MatrixBase<OtherDerived> const & other)
531  {
532  EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
533  call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
534  }
535 
536  template<typename RhsType, typename DstType>
537  EIGEN_DEVICE_FUNC
538  EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const {
539  if(!internal::is_same_dense(dst,rhs))
540  dst = rhs;
541  this->solveInPlace(dst);
542  }
543 
544  template<typename ProductType>
545  EIGEN_DEVICE_FUNC
546  EIGEN_STRONG_INLINE TriangularViewType& _assignProduct(const ProductType& prod, const Scalar& alpha, bool beta);
547 };
548 
549 /***************************************************************************
550 * Implementation of triangular evaluation/assignment
551 ***************************************************************************/
552 
553 #ifndef EIGEN_PARSED_BY_DOXYGEN
554 // FIXME should we keep that possibility
555 template<typename MatrixType, unsigned int Mode>
556 template<typename OtherDerived>
557 inline TriangularView<MatrixType, Mode>&
558 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const MatrixBase<OtherDerived>& other)
559 {
560  internal::call_assignment_no_alias(derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
561  return derived();
562 }
563 
564 // FIXME should we keep that possibility
565 template<typename MatrixType, unsigned int Mode>
566 template<typename OtherDerived>
567 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const MatrixBase<OtherDerived>& other)
568 {
569  internal::call_assignment_no_alias(derived(), other.template triangularView<Mode>());
570 }
571 
572 
573 
574 template<typename MatrixType, unsigned int Mode>
575 template<typename OtherDerived>
576 inline TriangularView<MatrixType, Mode>&
577 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const TriangularBase<OtherDerived>& other)
578 {
579  eigen_assert(Mode == int(OtherDerived::Mode));
580  internal::call_assignment(derived(), other.derived());
581  return derived();
582 }
583 
584 template<typename MatrixType, unsigned int Mode>
585 template<typename OtherDerived>
586 void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const TriangularBase<OtherDerived>& other)
587 {
588  eigen_assert(Mode == int(OtherDerived::Mode));
589  internal::call_assignment_no_alias(derived(), other.derived());
590 }
591 #endif
592 
593 /***************************************************************************
594 * Implementation of TriangularBase methods
595 ***************************************************************************/
596 
599 template<typename Derived>
600 template<typename DenseDerived>
602 {
603  evalToLazy(other.derived());
604 }
605 
606 /***************************************************************************
607 * Implementation of TriangularView methods
608 ***************************************************************************/
609 
610 /***************************************************************************
611 * Implementation of MatrixBase methods
612 ***************************************************************************/
613 
625 template<typename Derived>
626 template<unsigned int Mode>
627 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
629 {
630  return typename TriangularViewReturnType<Mode>::Type(derived());
631 }
632 
634 template<typename Derived>
635 template<unsigned int Mode>
636 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
638 {
639  return typename ConstTriangularViewReturnType<Mode>::Type(derived());
640 }
641 
647 template<typename Derived>
648 bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const
649 {
650  RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
651  for(Index j = 0; j < cols(); ++j)
652  {
653  Index maxi = numext::mini(j, rows()-1);
654  for(Index i = 0; i <= maxi; ++i)
655  {
656  RealScalar absValue = numext::abs(coeff(i,j));
657  if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
658  }
659  }
660  RealScalar threshold = maxAbsOnUpperPart * prec;
661  for(Index j = 0; j < cols(); ++j)
662  for(Index i = j+1; i < rows(); ++i)
663  if(numext::abs(coeff(i, j)) > threshold) return false;
664  return true;
665 }
666 
672 template<typename Derived>
673 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const
674 {
675  RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
676  for(Index j = 0; j < cols(); ++j)
677  for(Index i = j; i < rows(); ++i)
678  {
679  RealScalar absValue = numext::abs(coeff(i,j));
680  if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
681  }
682  RealScalar threshold = maxAbsOnLowerPart * prec;
683  for(Index j = 1; j < cols(); ++j)
684  {
685  Index maxi = numext::mini(j, rows()-1);
686  for(Index i = 0; i < maxi; ++i)
687  if(numext::abs(coeff(i, j)) > threshold) return false;
688  }
689  return true;
690 }
691 
692 
693 /***************************************************************************
694 ****************************************************************************
695 * Evaluators and Assignment of triangular expressions
696 ***************************************************************************
697 ***************************************************************************/
698 
699 namespace internal {
700 
701 
702 // TODO currently a triangular expression has the form TriangularView<.,.>
703 // in the future triangular-ness should be defined by the expression traits
704 // such that Transpose<TriangularView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
705 template<typename MatrixType, unsigned int Mode>
706 struct evaluator_traits<TriangularView<MatrixType,Mode> >
707 {
708  typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
709  typedef typename glue_shapes<typename evaluator_traits<MatrixType>::Shape, TriangularShape>::type Shape;
710 };
711 
712 template<typename MatrixType, unsigned int Mode>
713 struct unary_evaluator<TriangularView<MatrixType,Mode>, IndexBased>
714  : evaluator<typename internal::remove_all<MatrixType>::type>
715 {
716  typedef TriangularView<MatrixType,Mode> XprType;
717  typedef evaluator<typename internal::remove_all<MatrixType>::type> Base;
718  unary_evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {}
719 };
720 
721 // Additional assignment kinds:
722 struct Triangular2Triangular {};
723 struct Triangular2Dense {};
724 struct Dense2Triangular {};
725 
726 
727 template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop;
728 
729 
735 template<int UpLo, int Mode, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version = Specialized>
736 class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version>
737 {
738 protected:
739  typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base;
740  typedef typename Base::DstXprType DstXprType;
741  typedef typename Base::SrcXprType SrcXprType;
742  using Base::m_dst;
743  using Base::m_src;
744  using Base::m_functor;
745 public:
746 
747  typedef typename Base::DstEvaluatorType DstEvaluatorType;
748  typedef typename Base::SrcEvaluatorType SrcEvaluatorType;
749  typedef typename Base::Scalar Scalar;
750  typedef typename Base::AssignmentTraits AssignmentTraits;
751 
752 
753  EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr)
754  : Base(dst, src, func, dstExpr)
755  {}
756 
757 #ifdef EIGEN_INTERNAL_DEBUGGING
758  EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col)
759  {
760  eigen_internal_assert(row!=col);
761  Base::assignCoeff(row,col);
762  }
763 #else
764  using Base::assignCoeff;
765 #endif
766 
767  EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id)
768  {
769  if(Mode==UnitDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(1));
770  else if(Mode==ZeroDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(0));
771  else if(Mode==0) Base::assignCoeff(id,id);
772  }
773 
774  EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index row, Index col)
775  {
776  eigen_internal_assert(row!=col);
777  if(SetOpposite)
778  m_functor.assignCoeff(m_dst.coeffRef(row,col), Scalar(0));
779  }
780 };
781 
782 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType, typename Functor>
783 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
784 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src, const Functor &func)
785 {
786  typedef evaluator<DstXprType> DstEvaluatorType;
787  typedef evaluator<SrcXprType> SrcEvaluatorType;
788 
789  SrcEvaluatorType srcEvaluator(src);
790 
791  Index dstRows = src.rows();
792  Index dstCols = src.cols();
793  if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
794  dst.resize(dstRows, dstCols);
795  DstEvaluatorType dstEvaluator(dst);
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 * (DstEvaluatorType::CoeffReadCost+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 EIGEN_STRONG_INLINE
812 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src)
813 {
814  call_triangular_assignment_loop<Mode,SetOpposite>(dst, src, internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>());
815 }
816 
817 template<> struct AssignmentKind<TriangularShape,TriangularShape> { typedef Triangular2Triangular Kind; };
818 template<> struct AssignmentKind<DenseShape,TriangularShape> { typedef Triangular2Dense Kind; };
819 template<> struct AssignmentKind<TriangularShape,DenseShape> { typedef Dense2Triangular Kind; };
820 
821 
822 template< typename DstXprType, typename SrcXprType, typename Functor>
823 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Triangular>
824 {
825  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
826  {
827  eigen_assert(int(DstXprType::Mode) == int(SrcXprType::Mode));
828 
829  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
830  }
831 };
832 
833 template< typename DstXprType, typename SrcXprType, typename Functor>
834 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Dense>
835 {
836  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
837  {
838  call_triangular_assignment_loop<SrcXprType::Mode, (SrcXprType::Mode&SelfAdjoint)==0>(dst, src, func);
839  }
840 };
841 
842 template< typename DstXprType, typename SrcXprType, typename Functor>
843 struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular>
844 {
845  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
846  {
847  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
848  }
849 };
850 
851 
852 template<typename Kernel, unsigned int Mode, int UnrollCount, bool SetOpposite>
853 struct triangular_assignment_loop
854 {
855  // FIXME: this is not very clean, perhaps this information should be provided by the kernel?
856  typedef typename Kernel::DstEvaluatorType DstEvaluatorType;
857  typedef typename DstEvaluatorType::XprType DstXprType;
858 
859  enum {
860  col = (UnrollCount-1) / DstXprType::RowsAtCompileTime,
861  row = (UnrollCount-1) % DstXprType::RowsAtCompileTime
862  };
863 
864  typedef typename Kernel::Scalar Scalar;
865 
866  EIGEN_DEVICE_FUNC
867  static inline void run(Kernel &kernel)
868  {
869  triangular_assignment_loop<Kernel, Mode, UnrollCount-1, SetOpposite>::run(kernel);
870 
871  if(row==col)
872  kernel.assignDiagonalCoeff(row);
873  else if( ((Mode&Lower) && row>col) || ((Mode&Upper) && row<col) )
874  kernel.assignCoeff(row,col);
875  else if(SetOpposite)
876  kernel.assignOppositeCoeff(row,col);
877  }
878 };
879 
880 // prevent buggy user code from causing an infinite recursion
881 template<typename Kernel, unsigned int Mode, bool SetOpposite>
882 struct triangular_assignment_loop<Kernel, Mode, 0, SetOpposite>
883 {
884  EIGEN_DEVICE_FUNC
885  static inline void run(Kernel &) {}
886 };
887 
888 
889 
890 // TODO: experiment with a recursive assignment procedure splitting the current
891 // triangular part into one rectangular and two triangular parts.
892 
893 
894 template<typename Kernel, unsigned int Mode, bool SetOpposite>
895 struct triangular_assignment_loop<Kernel, Mode, Dynamic, SetOpposite>
896 {
897  typedef typename Kernel::Scalar Scalar;
898  EIGEN_DEVICE_FUNC
899  static inline void run(Kernel &kernel)
900  {
901  for(Index j = 0; j < kernel.cols(); ++j)
902  {
903  Index maxi = numext::mini(j, kernel.rows());
904  Index i = 0;
905  if (((Mode&Lower) && SetOpposite) || (Mode&Upper))
906  {
907  for(; i < maxi; ++i)
908  if(Mode&Upper) kernel.assignCoeff(i, j);
909  else kernel.assignOppositeCoeff(i, j);
910  }
911  else
912  i = maxi;
913 
914  if(i<kernel.rows()) // then i==j
915  kernel.assignDiagonalCoeff(i++);
916 
917  if (((Mode&Upper) && SetOpposite) || (Mode&Lower))
918  {
919  for(; i < kernel.rows(); ++i)
920  if(Mode&Lower) kernel.assignCoeff(i, j);
921  else kernel.assignOppositeCoeff(i, j);
922  }
923  }
924  }
925 };
926 
927 } // end namespace internal
928 
931 template<typename Derived>
932 template<typename DenseDerived>
934 {
935  other.derived().resize(this->rows(), this->cols());
936  internal::call_triangular_assignment_loop<Derived::Mode,(Derived::Mode&SelfAdjoint)==0 /* SetOpposite */>(other.derived(), derived().nestedExpression());
937 }
938 
939 namespace internal {
940 
941 // Triangular = Product
942 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
943 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
944 {
945  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
946  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,typename SrcXprType::Scalar> &)
947  {
948  Index dstRows = src.rows();
949  Index dstCols = src.cols();
950  if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
951  dst.resize(dstRows, dstCols);
952 
953  dst._assignProduct(src, 1, 0);
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::add_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::add_assign_op<Scalar,typename SrcXprType::Scalar> &)
963  {
964  dst._assignProduct(src, 1, 1);
965  }
966 };
967 
968 // Triangular -= Product
969 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
970 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::sub_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
971 {
972  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
973  static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,typename SrcXprType::Scalar> &)
974  {
975  dst._assignProduct(src, -1, 1);
976  }
977 };
978 
979 } // end namespace internal
980 
981 } // end namespace Eigen
982 
983 #endif // EIGEN_TRIANGULARMATRIX_H
const Product< TriangularViewType, OtherDerived > operator*(const MatrixBase< OtherDerived > &rhs) const
Definition: TriangularMatrix.h:453
Scalar determinant() const
Definition: TriangularMatrix.h:309
const int HugeCost
Definition: Constants.h:39
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:71
Index outerStride() const
Definition: TriangularMatrix.h:360
NestedExpression & nestedExpression()
Definition: TriangularMatrix.h:237
TriangularViewType & setConstant(const Scalar &value)
Definition: TriangularMatrix.h:393
Definition: Constants.h:216
TransposeReturnType transpose()
Definition: TriangularMatrix.h:254
Base class for triangular part in a matrix.
Definition: TriangularMatrix.h:27
const unsigned int DirectAccessBit
Definition: Constants.h:150
const AdjointReturnType adjoint() const
Definition: TriangularMatrix.h:248
const SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView() const
Definition: TriangularMatrix.h:299
const unsigned int LvalueBit
Definition: Constants.h:139
friend const Product< OtherDerived, TriangularViewType > operator*(const MatrixBase< OtherDerived > &lhs, const TriangularViewImpl &rhs)
Definition: TriangularMatrix.h:462
Namespace containing all symbols from the Eigen library.
Definition: Core:306
TriangularViewType & setOnes()
Definition: TriangularMatrix.h:400
Derived & derived()
Definition: EigenBase.h:45
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:38
TriangularViewType & operator*=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:383
Base class for all dense matrices, vectors, and arrays.
Definition: DenseBase.h:41
const unsigned int PacketAccessBit
Definition: Constants.h:89
const ConjugateReturnType conjugate() const
Definition: TriangularMatrix.h:242
Definition: EigenBase.h:29
bool isUpperTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:648
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:530
Derived & derived()
Definition: EigenBase.h:45
Definition: Constants.h:208
Scalar coeff(Index row, Index col) const
Definition: TriangularMatrix.h:406
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
void swap(TriangularBase< OtherDerived > &other)
Definition: TriangularMatrix.h:517
void evalTo(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:601
TriangularViewType & operator/=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:386
Expression of a selfadjoint matrix from a triangular part of a dense matrix.
Definition: SelfAdjointView.h:49
Definition: Constants.h:218
Definition: Constants.h:214
Definition: Constants.h:206
TriangularViewType & setZero()
Definition: TriangularMatrix.h:397
Definition: Constants.h:210
Definition: Constants.h:212
Definition: Eigen_Colamd.h:50
bool isLowerTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:673
Index rows() const
Definition: TriangularMatrix.h:226
Scalar & coeffRef(Index row, Index col)
Definition: TriangularMatrix.h:416
Definition: Constants.h:220
Index cols() const
Definition: TriangularMatrix.h:229
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:186
SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView()
Definition: TriangularMatrix.h:291
TriangularViewType & operator-=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:376
Definition: Constants.h:491
Index innerStride() const
Definition: TriangularMatrix.h:364
const int Dynamic
Definition: Constants.h:21
Pseudo expression representing a solving operation.
Definition: Solve.h:62
const NestedExpression & nestedExpression() const
Definition: TriangularMatrix.h:233
const ConstTransposeReturnType transpose() const
Definition: TriangularMatrix.h:264
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const unsigned int LinearAccessBit
Definition: Constants.h:125
void fill(const Scalar &value)
Definition: TriangularMatrix.h:390
void evalToLazy(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:933
TriangularViewType & operator+=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:369
Definition: TriangularMatrix.h:38