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