Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
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 #include "./InternalHeaderCheck.h"
15 
16 namespace Eigen {
17 
18 namespace internal {
19 
20 template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval;
21 
22 }
23 
29 template<typename Derived> class TriangularBase : public EigenBase<Derived>
30 {
31  public:
32 
33  enum {
34  Mode = internal::traits<Derived>::Mode,
35  RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
36  ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
37  MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
38  MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime,
39 
40  SizeAtCompileTime = (internal::size_of_xpr_at_compile_time<Derived>::ret),
45  MaxSizeAtCompileTime = internal::size_at_compile_time(internal::traits<Derived>::MaxRowsAtCompileTime,
46  internal::traits<Derived>::MaxColsAtCompileTime)
47 
48  };
49  typedef typename internal::traits<Derived>::Scalar Scalar;
50  typedef typename internal::traits<Derived>::StorageKind StorageKind;
51  typedef typename internal::traits<Derived>::StorageIndex StorageIndex;
52  typedef typename internal::traits<Derived>::FullMatrixType DenseMatrixType;
53  typedef DenseMatrixType DenseType;
54  typedef Derived const& Nested;
55 
56  EIGEN_DEVICE_FUNC
57  inline TriangularBase() { eigen_assert(!((int(Mode) & int(UnitDiag)) && (int(Mode) & int(ZeroDiag)))); }
58 
59  EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
60  inline Index rows() const EIGEN_NOEXCEPT { return derived().rows(); }
61  EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
62  inline Index cols() const EIGEN_NOEXCEPT { return derived().cols(); }
63  EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
64  inline Index outerStride() const EIGEN_NOEXCEPT { return derived().outerStride(); }
65  EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
66  inline Index innerStride() const EIGEN_NOEXCEPT { return derived().innerStride(); }
67 
68  // dummy resize function
69  EIGEN_DEVICE_FUNC
70  void resize(Index rows, Index cols)
71  {
72  EIGEN_UNUSED_VARIABLE(rows);
73  EIGEN_UNUSED_VARIABLE(cols);
74  eigen_assert(rows==this->rows() && cols==this->cols());
75  }
76 
77  EIGEN_DEVICE_FUNC
78  inline Scalar coeff(Index row, Index col) const { return derived().coeff(row,col); }
79  EIGEN_DEVICE_FUNC
80  inline Scalar& coeffRef(Index row, Index col) { return derived().coeffRef(row,col); }
81 
84  template<typename Other>
85  EIGEN_DEVICE_FUNC
86  EIGEN_STRONG_INLINE void copyCoeff(Index row, Index col, Other& other)
87  {
88  derived().coeffRef(row, col) = other.coeff(row, col);
89  }
90 
91  EIGEN_DEVICE_FUNC
92  inline Scalar operator()(Index row, Index col) const
93  {
94  check_coordinates(row, col);
95  return coeff(row,col);
96  }
97  EIGEN_DEVICE_FUNC
98  inline Scalar& operator()(Index row, Index col)
99  {
100  check_coordinates(row, col);
101  return coeffRef(row,col);
102  }
103 
104  #ifndef EIGEN_PARSED_BY_DOXYGEN
105  EIGEN_DEVICE_FUNC
106  inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
107  EIGEN_DEVICE_FUNC
108  inline Derived& derived() { return *static_cast<Derived*>(this); }
109  #endif // not EIGEN_PARSED_BY_DOXYGEN
110 
111  template<typename DenseDerived>
112  EIGEN_DEVICE_FUNC
113  void evalTo(MatrixBase<DenseDerived> &other) const;
114  template<typename DenseDerived>
115  EIGEN_DEVICE_FUNC
117 
118  EIGEN_DEVICE_FUNC
119  DenseMatrixType toDenseMatrix() const
120  {
121  DenseMatrixType res(rows(), cols());
122  evalToLazy(res);
123  return res;
124  }
125 
126  protected:
127 
128  void check_coordinates(Index row, Index col) const
129  {
130  EIGEN_ONLY_USED_FOR_DEBUG(row);
131  EIGEN_ONLY_USED_FOR_DEBUG(col);
132  eigen_assert(col>=0 && col<cols() && row>=0 && row<rows());
133  const int mode = int(Mode) & ~SelfAdjoint;
134  EIGEN_ONLY_USED_FOR_DEBUG(mode);
135  eigen_assert((mode==Upper && col>=row)
136  || (mode==Lower && col<=row)
137  || ((mode==StrictlyUpper || mode==UnitUpper) && col>row)
138  || ((mode==StrictlyLower || mode==UnitLower) && col<row));
139  }
140 
141  #ifdef EIGEN_INTERNAL_DEBUGGING
142  void check_coordinates_internal(Index row, Index col) const
143  {
144  check_coordinates(row, col);
145  }
146  #else
147  void check_coordinates_internal(Index , Index ) const {}
148  #endif
149 
150 };
151 
169 namespace internal {
170 template<typename MatrixType, unsigned int Mode_>
171 struct traits<TriangularView<MatrixType, Mode_> > : traits<MatrixType>
172 {
173  typedef typename ref_selector<MatrixType>::non_const_type MatrixTypeNested;
174  typedef std::remove_reference_t<MatrixTypeNested> MatrixTypeNestedNonRef;
175  typedef remove_all_t<MatrixTypeNested> MatrixTypeNestedCleaned;
176  typedef typename MatrixType::PlainObject FullMatrixType;
177  typedef MatrixType ExpressionType;
178  enum {
179  Mode = Mode_,
180  FlagsLvalueBit = is_lvalue<MatrixType>::value ? LvalueBit : 0,
181  Flags = (MatrixTypeNestedCleaned::Flags & (HereditaryBits | FlagsLvalueBit) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)))
182  };
183 };
184 }
185 
186 template<typename MatrixType_, unsigned int Mode_, typename StorageKind> class TriangularViewImpl;
187 
188 template<typename MatrixType_, unsigned int Mode_> class TriangularView
189  : public TriangularViewImpl<MatrixType_, Mode_, typename internal::traits<MatrixType_>::StorageKind >
190 {
191  public:
192 
193  typedef TriangularViewImpl<MatrixType_, Mode_, typename internal::traits<MatrixType_>::StorageKind > Base;
194  typedef typename internal::traits<TriangularView>::Scalar Scalar;
195  typedef MatrixType_ MatrixType;
196 
197  protected:
198  typedef typename internal::traits<TriangularView>::MatrixTypeNested MatrixTypeNested;
199  typedef typename internal::traits<TriangularView>::MatrixTypeNestedNonRef MatrixTypeNestedNonRef;
200 
201  typedef internal::remove_all_t<typename MatrixType::ConjugateReturnType> MatrixConjugateReturnType;
203 
204  public:
205 
206  typedef typename internal::traits<TriangularView>::StorageKind StorageKind;
207  typedef typename internal::traits<TriangularView>::MatrixTypeNestedCleaned NestedExpression;
208 
209  enum {
210  Mode = Mode_,
211  Flags = internal::traits<TriangularView>::Flags,
212  TransposeMode = (Mode & Upper ? Lower : 0)
213  | (Mode & Lower ? Upper : 0)
214  | (Mode & (UnitDiag))
215  | (Mode & (ZeroDiag)),
216  IsVectorAtCompileTime = false
217  };
218 
219  EIGEN_DEVICE_FUNC
220  explicit inline TriangularView(MatrixType& matrix) : m_matrix(matrix)
221  {}
222 
223  EIGEN_INHERIT_ASSIGNMENT_OPERATORS(TriangularView)
224 
225 
226  EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
227  inline Index rows() const EIGEN_NOEXCEPT { return m_matrix.rows(); }
229  EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
230  inline Index cols() const EIGEN_NOEXCEPT { 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 m_matrix; }
239 
240  typedef TriangularView<const MatrixConjugateReturnType,Mode> ConjugateReturnType;
242  EIGEN_DEVICE_FUNC
243  inline const ConjugateReturnType conjugate() const
244  { return ConjugateReturnType(m_matrix.conjugate()); }
245 
249  template<bool Cond>
250  EIGEN_DEVICE_FUNC
251  inline std::conditional_t<Cond,ConjugateReturnType,ConstTriangularView>
252  conjugateIf() const
253  {
254  typedef std::conditional_t<Cond,ConjugateReturnType,ConstTriangularView> ReturnType;
255  return ReturnType(m_matrix.template conjugateIf<Cond>());
256  }
257 
260  EIGEN_DEVICE_FUNC
261  inline const AdjointReturnType adjoint() const
262  { return AdjointReturnType(m_matrix.adjoint()); }
263 
266  template<class Dummy=int>
267  EIGEN_DEVICE_FUNC
268  inline TransposeReturnType transpose(std::enable_if_t<Eigen::internal::is_lvalue<MatrixType>::value, Dummy*> = nullptr)
269  {
270  typename MatrixType::TransposeReturnType tmp(m_matrix);
271  return TransposeReturnType(tmp);
272  }
273 
276  EIGEN_DEVICE_FUNC
277  inline const ConstTransposeReturnType transpose() const
278  {
279  return ConstTransposeReturnType(m_matrix.transpose());
280  }
281 
282  template<typename Other>
283  EIGEN_DEVICE_FUNC
284  inline const Solve<TriangularView, Other>
285  solve(const MatrixBase<Other>& other) const
286  { return Solve<TriangularView, Other>(*this, other.derived()); }
287 
288  // workaround MSVC ICE
289  #if EIGEN_COMP_MSVC
290  template<int Side, typename Other>
291  EIGEN_DEVICE_FUNC
292  inline const internal::triangular_solve_retval<Side,TriangularView, Other>
293  solve(const MatrixBase<Other>& other) const
294  { return Base::template solve<Side>(other); }
295  #else
296  using Base::solve;
297  #endif
298 
303  EIGEN_DEVICE_FUNC
305  {
306  EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
308  }
309 
311  EIGEN_DEVICE_FUNC
313  {
314  EIGEN_STATIC_ASSERT((Mode&(UnitDiag|ZeroDiag))==0,PROGRAMMING_ERROR);
316  }
317 
318 
321  EIGEN_DEVICE_FUNC
322  Scalar determinant() const
323  {
324  if (Mode & UnitDiag)
325  return 1;
326  else if (Mode & ZeroDiag)
327  return 0;
328  else
329  return m_matrix.diagonal().prod();
330  }
331 
332  protected:
333 
334  MatrixTypeNested m_matrix;
335 };
336 
346 template<typename MatrixType_, unsigned int Mode_> class TriangularViewImpl<MatrixType_,Mode_,Dense>
347  : public TriangularBase<TriangularView<MatrixType_, Mode_> >
348 {
349  public:
350 
352 
354  typedef typename internal::traits<TriangularViewType>::Scalar Scalar;
355 
356  typedef MatrixType_ MatrixType;
357  typedef typename MatrixType::PlainObject DenseMatrixType;
358  typedef DenseMatrixType PlainObject;
359 
360  public:
361  using Base::evalToLazy;
362  using Base::derived;
363 
364  typedef typename internal::traits<TriangularViewType>::StorageKind StorageKind;
365 
366  enum {
367  Mode = Mode_,
368  Flags = internal::traits<TriangularViewType>::Flags
369  };
370 
373  EIGEN_DEVICE_FUNC
374  inline Index outerStride() const { return derived().nestedExpression().outerStride(); }
377  EIGEN_DEVICE_FUNC
378  inline Index innerStride() const { return derived().nestedExpression().innerStride(); }
379 
381  template<typename Other>
382  EIGEN_DEVICE_FUNC
384  internal::call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op<Scalar,typename Other::Scalar>());
385  return derived();
386  }
388  template<typename Other>
389  EIGEN_DEVICE_FUNC
391  internal::call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op<Scalar,typename Other::Scalar>());
392  return derived();
393  }
394 
396  EIGEN_DEVICE_FUNC
397  TriangularViewType& operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() * other; }
399  EIGEN_DEVICE_FUNC
400  TriangularViewType& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = derived().nestedExpression() / other; }
401 
403  EIGEN_DEVICE_FUNC
404  void fill(const Scalar& value) { setConstant(value); }
406  EIGEN_DEVICE_FUNC
407  TriangularViewType& setConstant(const Scalar& value)
408  { return *this = MatrixType::Constant(derived().rows(), derived().cols(), value); }
410  EIGEN_DEVICE_FUNC
411  TriangularViewType& setZero() { return setConstant(Scalar(0)); }
413  EIGEN_DEVICE_FUNC
414  TriangularViewType& setOnes() { return setConstant(Scalar(1)); }
415 
419  EIGEN_DEVICE_FUNC
420  inline Scalar coeff(Index row, Index col) const
421  {
422  Base::check_coordinates_internal(row, col);
423  return derived().nestedExpression().coeff(row, col);
424  }
425 
429  EIGEN_DEVICE_FUNC
430  inline Scalar& coeffRef(Index row, Index col)
431  {
432  EIGEN_STATIC_ASSERT_LVALUE(TriangularViewType);
433  Base::check_coordinates_internal(row, col);
434  return derived().nestedExpression().coeffRef(row, col);
435  }
436 
438  template<typename OtherDerived>
439  EIGEN_DEVICE_FUNC
441 
443  template<typename OtherDerived>
444  EIGEN_DEVICE_FUNC
446 
447 #ifndef EIGEN_PARSED_BY_DOXYGEN
448  EIGEN_DEVICE_FUNC
449  TriangularViewType& operator=(const TriangularViewImpl& other)
450  { return *this = other.derived().nestedExpression(); }
451 
452  template<typename OtherDerived>
454  EIGEN_DEPRECATED EIGEN_DEVICE_FUNC
455  void lazyAssign(const TriangularBase<OtherDerived>& other);
456 
457  template<typename OtherDerived>
459  EIGEN_DEPRECATED EIGEN_DEVICE_FUNC
460  void lazyAssign(const MatrixBase<OtherDerived>& other);
461 #endif
462 
464  template<typename OtherDerived>
465  EIGEN_DEVICE_FUNC
468  {
469  return Product<TriangularViewType,OtherDerived>(derived(), rhs.derived());
470  }
471 
473  template<typename OtherDerived> friend
474  EIGEN_DEVICE_FUNC
476  operator*(const MatrixBase<OtherDerived>& lhs, const TriangularViewImpl& rhs)
477  {
478  return Product<OtherDerived,TriangularViewType>(lhs.derived(),rhs.derived());
479  }
480 
504  template<int Side, typename Other>
505  inline const internal::triangular_solve_retval<Side,TriangularViewType, Other>
506  solve(const MatrixBase<Other>& other) const;
507 
517  template<int Side, typename OtherDerived>
518  EIGEN_DEVICE_FUNC
519  void solveInPlace(const MatrixBase<OtherDerived>& other) const;
520 
521  template<typename OtherDerived>
522  EIGEN_DEVICE_FUNC
523  void solveInPlace(const MatrixBase<OtherDerived>& other) const
524  { return solveInPlace<OnTheLeft>(other); }
525 
527  template<typename OtherDerived>
528  EIGEN_DEVICE_FUNC
529 #ifdef EIGEN_PARSED_BY_DOXYGEN
531 #else
532  void swap(TriangularBase<OtherDerived> const & other)
533 #endif
534  {
535  EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
536  call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
537  }
538 
540  template<typename OtherDerived>
542  EIGEN_DEPRECATED EIGEN_DEVICE_FUNC
543  void swap(MatrixBase<OtherDerived> const & other)
544  {
545  EIGEN_STATIC_ASSERT_LVALUE(OtherDerived);
546  call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op<Scalar>());
547  }
548 
549  template<typename RhsType, typename DstType>
550  EIGEN_DEVICE_FUNC
551  EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const {
552  if(!internal::is_same_dense(dst,rhs))
553  dst = rhs;
554  this->solveInPlace(dst);
555  }
556 
557  template<typename ProductType>
558  EIGEN_DEVICE_FUNC
559  EIGEN_STRONG_INLINE TriangularViewType& _assignProduct(const ProductType& prod, const Scalar& alpha, bool beta);
560  protected:
561  EIGEN_DEFAULT_COPY_CONSTRUCTOR(TriangularViewImpl)
562  EIGEN_DEFAULT_EMPTY_CONSTRUCTOR_AND_DESTRUCTOR(TriangularViewImpl)
563 
564 };
565 
566 /***************************************************************************
567 * Implementation of triangular evaluation/assignment
568 ***************************************************************************/
569 
570 #ifndef EIGEN_PARSED_BY_DOXYGEN
571 // FIXME should we keep that possibility
572 template<typename MatrixType, unsigned int Mode>
573 template<typename OtherDerived>
574 EIGEN_DEVICE_FUNC inline TriangularView<MatrixType, Mode>&
575 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const MatrixBase<OtherDerived>& other)
576 {
577  internal::call_assignment_no_alias(derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
578  return derived();
579 }
580 
581 // FIXME should we keep that possibility
582 template<typename MatrixType, unsigned int Mode>
583 template<typename OtherDerived>
584 EIGEN_DEVICE_FUNC void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const MatrixBase<OtherDerived>& other)
585 {
586  internal::call_assignment_no_alias(derived(), other.template triangularView<Mode>());
587 }
588 
589 
590 
591 template<typename MatrixType, unsigned int Mode>
592 template<typename OtherDerived>
593 EIGEN_DEVICE_FUNC inline TriangularView<MatrixType, Mode>&
594 TriangularViewImpl<MatrixType, Mode, Dense>::operator=(const TriangularBase<OtherDerived>& other)
595 {
596  eigen_assert(Mode == int(OtherDerived::Mode));
597  internal::call_assignment(derived(), other.derived());
598  return derived();
599 }
600 
601 template<typename MatrixType, unsigned int Mode>
602 template<typename OtherDerived>
603 EIGEN_DEVICE_FUNC void TriangularViewImpl<MatrixType, Mode, Dense>::lazyAssign(const TriangularBase<OtherDerived>& other)
604 {
605  eigen_assert(Mode == int(OtherDerived::Mode));
606  internal::call_assignment_no_alias(derived(), other.derived());
607 }
608 #endif
609 
610 /***************************************************************************
611 * Implementation of TriangularBase methods
612 ***************************************************************************/
613 
616 template<typename Derived>
617 template<typename DenseDerived>
618 EIGEN_DEVICE_FUNC void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
619 {
620  evalToLazy(other.derived());
621 }
622 
623 /***************************************************************************
624 * Implementation of TriangularView methods
625 ***************************************************************************/
626 
627 /***************************************************************************
628 * Implementation of MatrixBase methods
629 ***************************************************************************/
630 
642 template<typename Derived>
643 template<unsigned int Mode>
644 EIGEN_DEVICE_FUNC
645 typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
647 {
648  return typename TriangularViewReturnType<Mode>::Type(derived());
649 }
650 
652 template<typename Derived>
653 template<unsigned int Mode>
654 EIGEN_DEVICE_FUNC
655 typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
657 {
658  return typename ConstTriangularViewReturnType<Mode>::Type(derived());
659 }
660 
666 template<typename Derived>
667 bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const
668 {
669  RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
670  for(Index j = 0; j < cols(); ++j)
671  {
672  Index maxi = numext::mini(j, rows()-1);
673  for(Index i = 0; i <= maxi; ++i)
674  {
675  RealScalar absValue = numext::abs(coeff(i,j));
676  if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
677  }
678  }
679  RealScalar threshold = maxAbsOnUpperPart * prec;
680  for(Index j = 0; j < cols(); ++j)
681  for(Index i = j+1; i < rows(); ++i)
682  if(numext::abs(coeff(i, j)) > threshold) return false;
683  return true;
684 }
685 
691 template<typename Derived>
692 bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const
693 {
694  RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
695  for(Index j = 0; j < cols(); ++j)
696  for(Index i = j; i < rows(); ++i)
697  {
698  RealScalar absValue = numext::abs(coeff(i,j));
699  if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
700  }
701  RealScalar threshold = maxAbsOnLowerPart * prec;
702  for(Index j = 1; j < cols(); ++j)
703  {
704  Index maxi = numext::mini(j, rows()-1);
705  for(Index i = 0; i < maxi; ++i)
706  if(numext::abs(coeff(i, j)) > threshold) return false;
707  }
708  return true;
709 }
710 
711 
712 /***************************************************************************
713 ****************************************************************************
714 * Evaluators and Assignment of triangular expressions
715 ***************************************************************************
716 ***************************************************************************/
717 
718 namespace internal {
719 
720 
721 // TODO currently a triangular expression has the form TriangularView<.,.>
722 // in the future triangular-ness should be defined by the expression traits
723 // such that Transpose<TriangularView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
724 template<typename MatrixType, unsigned int Mode>
725 struct evaluator_traits<TriangularView<MatrixType,Mode> >
726 {
727  typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
728  typedef typename glue_shapes<typename evaluator_traits<MatrixType>::Shape, TriangularShape>::type Shape;
729 };
730 
731 template<typename MatrixType, unsigned int Mode>
732 struct unary_evaluator<TriangularView<MatrixType,Mode>, IndexBased>
733  : evaluator<internal::remove_all_t<MatrixType>>
734 {
735  typedef TriangularView<MatrixType,Mode> XprType;
736  typedef evaluator<internal::remove_all_t<MatrixType>> Base;
737  EIGEN_DEVICE_FUNC
738  unary_evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {}
739 };
740 
741 // Additional assignment kinds:
742 struct Triangular2Triangular {};
743 struct Triangular2Dense {};
744 struct Dense2Triangular {};
745 
746 
747 template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop;
748 
749 
755 template<int UpLo, int Mode, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version = Specialized>
756 class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version>
757 {
758 protected:
759  typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base;
760  typedef typename Base::DstXprType DstXprType;
761  typedef typename Base::SrcXprType SrcXprType;
762  using Base::m_dst;
763  using Base::m_src;
764  using Base::m_functor;
765 public:
766 
767  typedef typename Base::DstEvaluatorType DstEvaluatorType;
768  typedef typename Base::SrcEvaluatorType SrcEvaluatorType;
769  typedef typename Base::Scalar Scalar;
770  typedef typename Base::AssignmentTraits AssignmentTraits;
771 
772 
773  EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr)
774  : Base(dst, src, func, dstExpr)
775  {}
776 
777 #ifdef EIGEN_INTERNAL_DEBUGGING
778  EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col)
779  {
780  eigen_internal_assert(row!=col);
781  Base::assignCoeff(row,col);
782  }
783 #else
784  using Base::assignCoeff;
785 #endif
786 
787  EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id)
788  {
789  if(Mode==UnitDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(1));
790  else if(Mode==ZeroDiag && SetOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(0));
791  else if(Mode==0) Base::assignCoeff(id,id);
792  }
793 
794  EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index row, Index col)
795  {
796  eigen_internal_assert(row!=col);
797  if(SetOpposite)
798  m_functor.assignCoeff(m_dst.coeffRef(row,col), Scalar(0));
799  }
800 };
801 
802 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType, typename Functor>
803 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
804 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src, const Functor &func)
805 {
806  typedef evaluator<DstXprType> DstEvaluatorType;
807  typedef evaluator<SrcXprType> SrcEvaluatorType;
808 
809  SrcEvaluatorType srcEvaluator(src);
810 
811  Index dstRows = src.rows();
812  Index dstCols = src.cols();
813  if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
814  dst.resize(dstRows, dstCols);
815  DstEvaluatorType dstEvaluator(dst);
816 
817  typedef triangular_dense_assignment_kernel< Mode&(Lower|Upper),Mode&(UnitDiag|ZeroDiag|SelfAdjoint),SetOpposite,
818  DstEvaluatorType,SrcEvaluatorType,Functor> Kernel;
819  Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived());
820 
821  enum {
822  unroll = DstXprType::SizeAtCompileTime != Dynamic
823  && SrcEvaluatorType::CoeffReadCost < HugeCost
824  && DstXprType::SizeAtCompileTime * (int(DstEvaluatorType::CoeffReadCost) + int(SrcEvaluatorType::CoeffReadCost)) / 2 <= EIGEN_UNROLLING_LIMIT
825  };
826 
827  triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, SetOpposite>::run(kernel);
828 }
829 
830 template<int Mode, bool SetOpposite, typename DstXprType, typename SrcXprType>
831 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
832 void call_triangular_assignment_loop(DstXprType& dst, const SrcXprType& src)
833 {
834  call_triangular_assignment_loop<Mode,SetOpposite>(dst, src, internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>());
835 }
836 
837 template<> struct AssignmentKind<TriangularShape,TriangularShape> { typedef Triangular2Triangular Kind; };
838 template<> struct AssignmentKind<DenseShape,TriangularShape> { typedef Triangular2Dense Kind; };
839 template<> struct AssignmentKind<TriangularShape,DenseShape> { typedef Dense2Triangular Kind; };
840 
841 
842 template< typename DstXprType, typename SrcXprType, typename Functor>
843 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Triangular>
844 {
845  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
846  {
847  eigen_assert(int(DstXprType::Mode) == int(SrcXprType::Mode));
848 
849  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
850  }
851 };
852 
853 template< typename DstXprType, typename SrcXprType, typename Functor>
854 struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Dense>
855 {
856  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
857  {
858  call_triangular_assignment_loop<SrcXprType::Mode, (int(SrcXprType::Mode) & int(SelfAdjoint)) == 0>(dst, src, func);
859  }
860 };
861 
862 template< typename DstXprType, typename SrcXprType, typename Functor>
863 struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular>
864 {
865  EIGEN_DEVICE_FUNC static void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
866  {
867  call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func);
868  }
869 };
870 
871 
872 template<typename Kernel, unsigned int Mode, int UnrollCount, bool SetOpposite>
873 struct triangular_assignment_loop
874 {
875  // FIXME: this is not very clean, perhaps this information should be provided by the kernel?
876  typedef typename Kernel::DstEvaluatorType DstEvaluatorType;
877  typedef typename DstEvaluatorType::XprType DstXprType;
878 
879  enum {
880  col = (UnrollCount-1) / DstXprType::RowsAtCompileTime,
881  row = (UnrollCount-1) % DstXprType::RowsAtCompileTime
882  };
883 
884  typedef typename Kernel::Scalar Scalar;
885 
886  EIGEN_DEVICE_FUNC
887  static inline void run(Kernel &kernel)
888  {
889  triangular_assignment_loop<Kernel, Mode, UnrollCount-1, SetOpposite>::run(kernel);
890 
891  if(row==col)
892  kernel.assignDiagonalCoeff(row);
893  else if( ((Mode&Lower) && row>col) || ((Mode&Upper) && row<col) )
894  kernel.assignCoeff(row,col);
895  else if(SetOpposite)
896  kernel.assignOppositeCoeff(row,col);
897  }
898 };
899 
900 // prevent buggy user code from causing an infinite recursion
901 template<typename Kernel, unsigned int Mode, bool SetOpposite>
902 struct triangular_assignment_loop<Kernel, Mode, 0, SetOpposite>
903 {
904  EIGEN_DEVICE_FUNC
905  static inline void run(Kernel &) {}
906 };
907 
908 
909 
910 // TODO: experiment with a recursive assignment procedure splitting the current
911 // triangular part into one rectangular and two triangular parts.
912 
913 
914 template<typename Kernel, unsigned int Mode, bool SetOpposite>
915 struct triangular_assignment_loop<Kernel, Mode, Dynamic, SetOpposite>
916 {
917  typedef typename Kernel::Scalar Scalar;
918  EIGEN_DEVICE_FUNC
919  static inline void run(Kernel &kernel)
920  {
921  for(Index j = 0; j < kernel.cols(); ++j)
922  {
923  Index maxi = numext::mini(j, kernel.rows());
924  Index i = 0;
925  if (((Mode&Lower) && SetOpposite) || (Mode&Upper))
926  {
927  for(; i < maxi; ++i)
928  if(Mode&Upper) kernel.assignCoeff(i, j);
929  else kernel.assignOppositeCoeff(i, j);
930  }
931  else
932  i = maxi;
933 
934  if(i<kernel.rows()) // then i==j
935  kernel.assignDiagonalCoeff(i++);
936 
937  if (((Mode&Upper) && SetOpposite) || (Mode&Lower))
938  {
939  for(; i < kernel.rows(); ++i)
940  if(Mode&Lower) kernel.assignCoeff(i, j);
941  else kernel.assignOppositeCoeff(i, j);
942  }
943  }
944  }
945 };
946 
947 } // end namespace internal
948 
951 template<typename Derived>
952 template<typename DenseDerived>
954 {
955  other.derived().resize(this->rows(), this->cols());
956  internal::call_triangular_assignment_loop<Derived::Mode, (int(Derived::Mode) & int(SelfAdjoint)) == 0 /* SetOpposite */>(other.derived(), derived().nestedExpression());
957 }
958 
959 namespace internal {
960 
961 // Triangular = Product
962 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
963 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
964 {
965  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
966  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,typename SrcXprType::Scalar> &)
967  {
968  Index dstRows = src.rows();
969  Index dstCols = src.cols();
970  if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
971  dst.resize(dstRows, dstCols);
972 
973  dst._assignProduct(src, Scalar(1), false);
974  }
975 };
976 
977 // Triangular += Product
978 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
979 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::add_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
980 {
981  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
982  static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar,typename SrcXprType::Scalar> &)
983  {
984  dst._assignProduct(src, Scalar(1), true);
985  }
986 };
987 
988 // Triangular -= Product
989 template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
990 struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::sub_assign_op<Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, Dense2Triangular>
991 {
992  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
993  static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,typename SrcXprType::Scalar> &)
994  {
995  dst._assignProduct(src, Scalar(-1), true);
996  }
997 };
998 
999 } // end namespace internal
1000 
1001 } // end namespace Eigen
1002 
1003 #endif // EIGEN_TRIANGULARMATRIX_H
Base class for all dense matrices, vectors, and arrays.
Definition: DenseBase.h:42
Derived & derived()
Definition: EigenBase.h:48
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
bool isLowerTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:692
bool isUpperTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Definition: TriangularMatrix.h:667
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:77
Expression of a selfadjoint matrix from a triangular part of a dense matrix.
Definition: SelfAdjointView.h:53
Pseudo expression representing a solving operation.
Definition: Solve.h:65
Expression of the transpose of a matrix.
Definition: Transpose.h:56
Base class for triangular part in a matrix.
Definition: TriangularMatrix.h:30
void copyCoeff(Index row, Index col, Other &other)
Definition: TriangularMatrix.h:86
@ SizeAtCompileTime
Definition: TriangularMatrix.h:40
void evalTo(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:618
void evalToLazy(MatrixBase< DenseDerived > &other) const
Definition: TriangularMatrix.h:953
TriangularViewType & setZero()
Definition: TriangularMatrix.h:411
Scalar coeff(Index row, Index col) const
Definition: TriangularMatrix.h:420
TriangularViewType & operator/=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:400
void fill(const Scalar &value)
Definition: TriangularMatrix.h:404
TriangularViewType & operator-=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:390
Index outerStride() const
Definition: TriangularMatrix.h:374
const Product< TriangularViewType, OtherDerived > operator*(const MatrixBase< OtherDerived > &rhs) const
Definition: TriangularMatrix.h:467
EIGEN_DEPRECATED void swap(MatrixBase< OtherDerived > const &other)
Definition: TriangularMatrix.h:543
TriangularViewType & operator=(const MatrixBase< OtherDerived > &other)
Index innerStride() const
Definition: TriangularMatrix.h:378
TriangularViewType & setOnes()
Definition: TriangularMatrix.h:414
TriangularViewType & operator*=(const typename internal::traits< MatrixType >::Scalar &other)
Definition: TriangularMatrix.h:397
TriangularViewType & operator=(const TriangularBase< OtherDerived > &other)
TriangularViewType & setConstant(const Scalar &value)
Definition: TriangularMatrix.h:407
const internal::triangular_solve_retval< Side, TriangularViewType, Other > solve(const MatrixBase< Other > &other) const
friend const Product< OtherDerived, TriangularViewType > operator*(const MatrixBase< OtherDerived > &lhs, const TriangularViewImpl &rhs)
Definition: TriangularMatrix.h:476
void solveInPlace(const MatrixBase< OtherDerived > &other) const
TriangularViewType & operator+=(const DenseBase< Other > &other)
Definition: TriangularMatrix.h:383
Scalar & coeffRef(Index row, Index col)
Definition: TriangularMatrix.h:430
void swap(TriangularBase< OtherDerived > &other)
Definition: TriangularMatrix.h:530
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:190
const ConjugateReturnType conjugate() const
Definition: TriangularMatrix.h:243
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: TriangularMatrix.h:227
Scalar determinant() const
Definition: TriangularMatrix.h:322
TransposeReturnType transpose(std::enable_if_t< Eigen::internal::is_lvalue< MatrixType >::value, Dummy * >=nullptr)
Definition: TriangularMatrix.h:268
NestedExpression & nestedExpression()
Definition: TriangularMatrix.h:238
const SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView() const
Definition: TriangularMatrix.h:312
const ConstTransposeReturnType transpose() const
Definition: TriangularMatrix.h:277
std::conditional_t< Cond, ConjugateReturnType, ConstTriangularView > conjugateIf() const
Definition: TriangularMatrix.h:252
const AdjointReturnType adjoint() const
Definition: TriangularMatrix.h:261
const NestedExpression & nestedExpression() const
Definition: TriangularMatrix.h:234
SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView()
Definition: TriangularMatrix.h:304
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: TriangularMatrix.h:230
@ StrictlyLower
Definition: Constants.h:223
@ UnitDiag
Definition: Constants.h:215
@ StrictlyUpper
Definition: Constants.h:225
@ UnitLower
Definition: Constants.h:219
@ ZeroDiag
Definition: Constants.h:217
@ SelfAdjoint
Definition: Constants.h:227
@ UnitUpper
Definition: Constants.h:221
@ Lower
Definition: Constants.h:211
@ Upper
Definition: Constants.h:213
const unsigned int PacketAccessBit
Definition: Constants.h:96
const unsigned int LinearAccessBit
Definition: Constants.h:132
const unsigned int DirectAccessBit
Definition: Constants.h:157
const unsigned int LvalueBit
Definition: Constants.h:146
Namespace containing all symbols from the Eigen library.
Definition: Core:139
const int HugeCost
Definition: Constants.h:46
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const int Dynamic
Definition: Constants.h:24
Definition: Constants.h:509
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41