TriangularMatrix.h
Go to the documentation of this file.
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  };
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
69  {
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
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 
127  {
130  eigen_assert(col>=0 && col<cols() && row>=0 && row<rows());
131  const int mode = int(Mode) & ~SelfAdjoint;
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
146  #endif
147 
148 };
149 
167 namespace internal {
168 template<typename MatrixType, unsigned int _Mode>
169 struct traits<TriangularView<MatrixType, _Mode> > : traits<MatrixType>
170 {
174  typedef typename MatrixType::PlainObject FullMatrixType;
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 
193  typedef _MatrixType MatrixType;
194 
195  protected:
198 
200 
201  public:
202 
205 
206  enum {
207  Mode = _Mode,
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=;
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 
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
254  inline TransposeReturnType transpose()
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
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 
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 
351 
352  enum {
353  Mode = _Mode,
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
369  TriangularViewType& operator+=(const DenseBase<Other>& other) {
371  return derived();
372  }
374  template<typename Other>
375  EIGEN_DEVICE_FUNC
376  TriangularViewType& operator-=(const DenseBase<Other>& other) {
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
454  {
455  return Product<TriangularViewType,OtherDerived>(derived(), rhs.derived());
456  }
457 
459  template<typename OtherDerived> friend
460  EIGEN_DEVICE_FUNC
463  {
464  return Product<OtherDerived,TriangularViewType>(lhs.derived(),rhs.derived());
465  }
466 
490  template<int Side, typename Other>
491  EIGEN_DEVICE_FUNC
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);
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>
559 {
561  return derived();
562 }
563 
564 // FIXME should we keep that possibility
565 template<typename MatrixType, unsigned int Mode>
566 template<typename OtherDerived>
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>
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>
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>
629 {
630  return typename TriangularViewReturnType<Mode>::Type(derived());
631 }
632 
634 template<typename Derived>
635 template<unsigned int Mode>
638 {
639  return typename ConstTriangularViewReturnType<Mode>::Type(derived());
640 }
641 
647 template<typename Derived>
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>
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 {
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 {
718  unary_evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {}
719 };
720 
721 // Additional assignment kinds:
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:
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 
749  typedef typename Base::Scalar Scalar;
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 
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 
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>
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  {
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 {
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 {
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 {
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
EIGEN_DEVICE_FUNC const ConstTransposeReturnType transpose() const
EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index row, Index col)
SCALAR Scalar
Definition: bench_gemm.cpp:33
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_assignment_no_alias(Dst &dst, const Src &src, const Func &func)
#define EIGEN_STRONG_INLINE
Definition: Macros.h:494
remove_reference< MatrixTypeNested >::type MatrixTypeNestedNonRef
EIGEN_DEVICE_FUNC Scalar & coeffRef(Index row, Index col)
internal::remove_all< typename MatrixType::ConjugateReturnType >::type MatrixConjugateReturnType
EIGEN_DEVICE_FUNC void solveInPlace(const MatrixBase< OtherDerived > &other) const
const int HugeCost
Definition: Constants.h:39
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:71
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index cols() const
Definition: Product.h:101
EIGEN_DEVICE_FUNC Scalar & coeffRef(Index row, Index col)
return int(ret)+1
EIGEN_DEVICE_FUNC void swap(TriangularBase< OtherDerived > const &other)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index rows() const
Definition: Product.h:100
void check_coordinates_internal(Index, Index) const
bool is_same_dense(const T1 &mat1, const T2 &mat2, typename enable_if< has_direct_access< T1 >::ret &&has_direct_access< T2 >::ret, T1 >::type *=0)
Definition: XprHelper.h:661
TriangularView< typename MatrixType::TransposeReturnType, TransposeMode > TransposeReturnType
EIGEN_DEVICE_FUNC void evalToLazy(MatrixBase< DenseDerived > &other) const
TriangularView< _MatrixType, _Mode > TriangularViewType
Base class for triangular part in a matrix.
const unsigned int DirectAccessBit
Definition: Constants.h:150
generic_dense_assignment_kernel< DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version > Base
const unsigned int LvalueBit
Definition: Constants.h:139
EIGEN_DEVICE_FUNC Index rows() const
EIGEN_DEVICE_FUNC TriangularViewType & operator+=(const DenseBase< Other > &other)
bool isLowerTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
internal::traits< Derived >::FullMatrixType DenseMatrixType
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
EIGEN_DEVICE_FUNC TriangularViewType & setOnes()
MatrixXf MatrixType
TriangularView< const typename MatrixType::AdjointReturnType, TransposeMode > AdjointReturnType
glue_shapes< typename evaluator_traits< MatrixType >::Shape, TriangularShape >::type Shape
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
Definition: StaticAssert.h:124
EIGEN_DEVICE_FUNC void resize(Index newSize)
Definition: DenseBase.h:241
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
EIGEN_DEVICE_FUNC TriangularViewType & operator/=(const typename internal::traits< MatrixType >::Scalar &other)
ref_selector< MatrixType >::non_const_type MatrixTypeNested
EIGEN_DEVICE_FUNC Index innerStride() const
MatrixTypeNested m_matrix
EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id)
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:38
EIGEN_DEVICE_FUNC TransposeReturnType transpose()
Base class for all dense matrices, vectors, and arrays.
Definition: DenseBase.h:41
static EIGEN_DEVICE_FUNC void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
EIGEN_DEVICE_FUNC Index outerStride() const
const unsigned int PacketAccessBit
Definition: Constants.h:89
#define EIGEN_STATIC_ASSERT_LVALUE(Derived)
Definition: StaticAssert.h:199
static EIGEN_DEVICE_FUNC void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
DenseMatrixType DenseType
friend EIGEN_DEVICE_FUNC const Product< OtherDerived, TriangularViewType > operator*(const MatrixBase< OtherDerived > &lhs, const TriangularViewImpl &rhs)
EIGEN_DEVICE_FUNC const NestedExpression & nestedExpression() const
void resize(Index rows, Index cols)
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
EIGEN_DEVICE_FUNC void evalTo(MatrixBase< DenseDerived > &other) const
internal::traits< TriangularView >::Scalar Scalar
EIGEN_DEVICE_FUNC Derived & const_cast_derived() const
Definition: EigenBase.h:51
TriangularView< const typename MatrixType::ConstTransposeReturnType, TransposeMode > ConstTransposeReturnType
Derived const & Nested
const unsigned int HereditaryBits
Definition: Constants.h:190
EIGEN_DEVICE_FUNC TriangularViewType & setConstant(const Scalar &value)
EIGEN_DEVICE_FUNC const Product< TriangularViewType, OtherDerived > operator*(const MatrixBase< OtherDerived > &rhs) const
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE internal::enable_if< NumTraits< T >::IsSigned||NumTraits< T >::IsComplex, typename NumTraits< T >::Real >::type abs(const T &x)
EIGEN_DEVICE_FUNC TriangularViewType & operator*=(const typename internal::traits< MatrixType >::Scalar &other)
EIGEN_DEVICE_FUNC const AdjointReturnType adjoint() const
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void copyCoeff(Index row, Index col, Other &other)
EIGEN_DEVICE_FUNC TriangularBase()
EIGEN_DEVICE_FUNC const Solve< TriangularView, Other > solve(const MatrixBase< Other > &other) const
m row(1)
EIGEN_DEVICE_FUNC const SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView() const
TriangularView< const MatrixConjugateReturnType, Mode > ConjugateReturnType
EIGEN_DEVICE_FUNC void swap(MatrixBase< OtherDerived > const &other)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
storage_kind_to_evaluator_kind< typename MatrixType::StorageKind >::Kind Kind
EIGEN_DEVICE_FUNC Scalar operator()(Index row, Index col) const
#define eigen_assert(x)
Definition: Macros.h:579
TriangularViewImpl< _MatrixType, _Mode, typename internal::traits< _MatrixType >::StorageKind > Base
EIGEN_DEVICE_FUNC const ConjugateReturnType conjugate() const
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_triangular_assignment_loop(DstXprType &dst, const SrcXprType &src, const Functor &func)
EIGEN_DEVICE_FUNC TriangularViewReturnType< Mode >::Type triangularView()
RealScalar alpha
internal::traits< TriangularView >::MatrixTypeNested MatrixTypeNested
EIGEN_DEVICE_FUNC DenseMatrixType toDenseMatrix() const
Expression of a selfadjoint matrix from a triangular part of a dense matrix.
EIGEN_DEVICE_FUNC const Derived & derived() const
EIGEN_DEVICE_FUNC Index cols() const
EIGEN_DEVICE_FUNC Derived & derived()
evaluator< typename internal::remove_all< MatrixType >::type > Base
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T mini(const T &x, const T &y)
internal::traits< TriangularView >::MatrixTypeNestedCleaned NestedExpression
EIGEN_DEVICE_FUNC NestedExpression & nestedExpression()
EIGEN_DEVICE_FUNC TriangularViewType & operator-=(const DenseBase< Other > &other)
EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType &dstExpr)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const
EIGEN_DEVICE_FUNC Scalar & operator()(Index row, Index col)
EIGEN_DEVICE_FUNC Index cols() const
v setConstant(3, 5)
int func(const int &a)
Definition: testDSF.cpp:225
EIGEN_DEVICE_FUNC Index rows() const
NumTraits< Scalar >::Real RealScalar
Definition: DenseBase.h:73
static EIGEN_DEVICE_FUNC void run(Kernel &kernel)
DenseIndex ret
Definition: level1_impl.h:59
EIGEN_DEVICE_FUNC Scalar determinant() const
bool isUpperTriangular(const RealScalar &prec=NumTraits< Scalar >::dummy_precision()) const
Expression of a triangular part in a matrix.
internal::traits< Derived >::StorageKind StorageKind
m col(1)
EIGEN_DEVICE_FUNC Scalar coeff(Index row, Index col) const
internal::traits< TriangularViewType >::StorageKind StorageKind
internal::conditional< NumTraits< Scalar >::IsComplex, const CwiseUnaryOp< internal::scalar_conjugate_op< Scalar >, const Derived >, const Derived & >::type ConjugateReturnType
const int Dynamic
Definition: Constants.h:21
Pseudo expression representing a solving operation.
Definition: Solve.h:62
static EIGEN_DEVICE_FUNC void run(DstXprType &dst, const SrcXprType &src, const Functor &func)
#define eigen_internal_assert(x)
Definition: Macros.h:585
internal::traits< Derived >::Scalar Scalar
EIGEN_DEVICE_FUNC Scalar coeff(Index row, Index col) const
EIGEN_DEVICE_FUNC TriangularViewType & operator=(const TriangularViewImpl &other)
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
EIGEN_DEVICE_FUNC void fill(const Scalar &value)
TriangularView & operator=(const TriangularView &other)
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
EIGEN_DEVICE_FUNC TriangularView(MatrixType &matrix)
const unsigned int LinearAccessBit
Definition: Constants.h:125
EIGEN_DEVICE_FUNC SelfAdjointView< MatrixTypeNestedNonRef, Mode > selfadjointView()
internal::traits< TriangularView >::StorageKind StorageKind
void check_coordinates(Index row, Index col) const
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_assignment(Dst &dst, const Src &src)
std::ptrdiff_t j
#define EIGEN_UNROLLING_LIMIT
Definition: Settings.h:24
#define EIGEN_UNUSED_VARIABLE(var)
Definition: Macros.h:618
EIGEN_DEVICE_FUNC TriangularViewType & setZero()
Definition: pytypes.h:897
void swap(scoped_array< T > &a, scoped_array< T > &b)
Definition: Memory.h:602
internal::traits< Derived >::StorageIndex StorageIndex
#define EIGEN_ONLY_USED_FOR_DEBUG(x)
Definition: Macros.h:591
internal::traits< TriangularView >::MatrixTypeNestedNonRef MatrixTypeNestedNonRef
const Product< Lhs, Rhs > prod(const Lhs &lhs, const Rhs &rhs)
Definition: evaluators.cpp:8
internal::traits< TriangularViewType >::Scalar Scalar


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:51:15