SolveTriangular.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-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_SOLVETRIANGULAR_H
11 #define EIGEN_SOLVETRIANGULAR_H
12 
13 namespace Eigen {
14 
15 namespace internal {
16 
17 // Forward declarations:
18 // The following two routines are implemented in the products/TriangularSolver*.h files
19 template<typename LhsScalar, typename RhsScalar, typename Index, int Side, int Mode, bool Conjugate, int StorageOrder>
21 
22 template <typename Scalar, typename Index, int Side, int Mode, bool Conjugate, int TriStorageOrder, int OtherStorageOrder, int OtherInnerStride>
24 
25 // small helper struct extracting some traits on the underlying solver operation
26 template<typename Lhs, typename Rhs, int Side>
28 {
29  private:
30  enum {
31  RhsIsVectorAtCompileTime = (Side==OnTheLeft ? Rhs::ColsAtCompileTime : Rhs::RowsAtCompileTime)==1
32  };
33  public:
34  enum {
35  Unrolling = (RhsIsVectorAtCompileTime && Rhs::SizeAtCompileTime != Dynamic && Rhs::SizeAtCompileTime <= 8)
38  };
39 };
40 
41 template<typename Lhs, typename Rhs,
42  int Side, // can be OnTheLeft/OnTheRight
43  int Mode, // can be Upper/Lower | UnitDiag
46  >
48 
49 template<typename Lhs, typename Rhs, int Side, int Mode>
51 {
52  typedef typename Lhs::Scalar LhsScalar;
53  typedef typename Rhs::Scalar RhsScalar;
57  static EIGEN_DEVICE_FUNC void run(const Lhs& lhs, Rhs& rhs)
58  {
59  ActualLhsType actualLhs = LhsProductTraits::extract(lhs);
60 
61  // FIXME find a way to allow an inner stride if packet_traits<Scalar>::size==1
62 
63  bool useRhsDirectly = Rhs::InnerStrideAtCompileTime==1 || rhs.innerStride()==1;
64 
65  ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhs,rhs.size(),
66  (useRhsDirectly ? rhs.data() : 0));
67 
68  if(!useRhsDirectly)
69  MappedRhs(actualRhs,rhs.size()) = rhs;
70 
71  triangular_solve_vector<LhsScalar, RhsScalar, Index, Side, Mode, LhsProductTraits::NeedToConjugate,
72  (int(Lhs::Flags) & RowMajorBit) ? RowMajor : ColMajor>
73  ::run(actualLhs.cols(), actualLhs.data(), actualLhs.outerStride(), actualRhs);
74 
75  if(!useRhsDirectly)
76  rhs = MappedRhs(actualRhs, rhs.size());
77  }
78 };
79 
80 // the rhs is a matrix
81 template<typename Lhs, typename Rhs, int Side, int Mode>
83 {
84  typedef typename Rhs::Scalar Scalar;
87 
88  static EIGEN_DEVICE_FUNC void run(const Lhs& lhs, Rhs& rhs)
89  {
90  typename internal::add_const_on_value_type<ActualLhsType>::type actualLhs = LhsProductTraits::extract(lhs);
91 
92  const Index size = lhs.rows();
93  const Index othersize = Side==OnTheLeft? rhs.cols() : rhs.rows();
94 
95  typedef internal::gemm_blocking_space<(Rhs::Flags&RowMajorBit) ? RowMajor : ColMajor,Scalar,Scalar,
96  Rhs::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime, Lhs::MaxRowsAtCompileTime,4> BlockingType;
97 
98  BlockingType blocking(rhs.rows(), rhs.cols(), size, 1, false);
99 
100  triangular_solve_matrix<Scalar,Index,Side,Mode,LhsProductTraits::NeedToConjugate,(int(Lhs::Flags) & RowMajorBit) ? RowMajor : ColMajor,
101  (Rhs::Flags&RowMajorBit) ? RowMajor : ColMajor, Rhs::InnerStrideAtCompileTime>
102  ::run(size, othersize, &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &rhs.coeffRef(0,0), rhs.innerStride(), rhs.outerStride(), blocking);
103  }
104 };
105 
106 /***************************************************************************
107 * meta-unrolling implementation
108 ***************************************************************************/
109 
110 template<typename Lhs, typename Rhs, int Mode, int LoopIndex, int Size,
111  bool Stop = LoopIndex==Size>
113 
114 template<typename Lhs, typename Rhs, int Mode, int LoopIndex, int Size>
115 struct triangular_solver_unroller<Lhs,Rhs,Mode,LoopIndex,Size,false> {
116  enum {
117  IsLower = ((Mode&Lower)==Lower),
118  DiagIndex = IsLower ? LoopIndex : Size - LoopIndex - 1,
119  StartIndex = IsLower ? 0 : DiagIndex+1
120  };
121  static EIGEN_DEVICE_FUNC void run(const Lhs& lhs, Rhs& rhs)
122  {
123  if (LoopIndex>0)
124  rhs.coeffRef(DiagIndex) -= lhs.row(DiagIndex).template segment<LoopIndex>(StartIndex).transpose()
125  .cwiseProduct(rhs.template segment<LoopIndex>(StartIndex)).sum();
126 
127  if(!(Mode & UnitDiag))
128  rhs.coeffRef(DiagIndex) /= lhs.coeff(DiagIndex,DiagIndex);
129 
131  }
132 };
133 
134 template<typename Lhs, typename Rhs, int Mode, int LoopIndex, int Size>
135 struct triangular_solver_unroller<Lhs,Rhs,Mode,LoopIndex,Size,true> {
136  static EIGEN_DEVICE_FUNC void run(const Lhs&, Rhs&) {}
137 };
138 
139 template<typename Lhs, typename Rhs, int Mode>
141  static EIGEN_DEVICE_FUNC void run(const Lhs& lhs, Rhs& rhs)
143 };
144 
145 template<typename Lhs, typename Rhs, int Mode>
147  static EIGEN_DEVICE_FUNC void run(const Lhs& lhs, Rhs& rhs)
148  {
149  Transpose<const Lhs> trLhs(lhs);
150  Transpose<Rhs> trRhs(rhs);
151 
153  ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag),
154  0,Rhs::SizeAtCompileTime>::run(trLhs,trRhs);
155  }
156 };
157 
158 } // end namespace internal
159 
160 /***************************************************************************
161 * TriangularView methods
162 ***************************************************************************/
163 
164 #ifndef EIGEN_PARSED_BY_DOXYGEN
165 template<typename MatrixType, unsigned int Mode>
166 template<int Side, typename OtherDerived>
168 {
169  OtherDerived& other = _other.const_cast_derived();
170  eigen_assert( derived().cols() == derived().rows() && ((Side==OnTheLeft && derived().cols() == other.rows()) || (Side==OnTheRight && derived().cols() == other.cols())) );
171  eigen_assert((!(int(Mode) & int(ZeroDiag))) && bool(int(Mode) & (int(Upper) | int(Lower))));
172  // If solving for a 0x0 matrix, nothing to do, simply return.
173  if (derived().cols() == 0)
174  return;
175 
176  enum { copy = (internal::traits<OtherDerived>::Flags & RowMajorBit) && OtherDerived::IsVectorAtCompileTime && OtherDerived::SizeAtCompileTime!=1};
177  typedef typename internal::conditional<copy,
179  OtherCopy otherCopy(other);
180 
182  Side, Mode>::run(derived().nestedExpression(), otherCopy);
183 
184  if (copy)
185  other = otherCopy;
186 }
187 
188 template<typename Derived, unsigned int Mode>
189 template<int Side, typename Other>
192 {
193  return internal::triangular_solve_retval<Side,TriangularViewType,Other>(derived(), other.derived());
194 }
195 #endif
196 
197 namespace internal {
198 
199 
200 template<int Side, typename TriangularType, typename Rhs>
201 struct traits<triangular_solve_retval<Side, TriangularType, Rhs> >
202 {
204 };
205 
206 template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval
207  : public ReturnByValue<triangular_solve_retval<Side, TriangularType, Rhs> >
208 {
211 
212  triangular_solve_retval(const TriangularType& tri, const Rhs& rhs)
213  : m_triangularMatrix(tri), m_rhs(rhs)
214  {}
215 
216  inline EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_rhs.rows(); }
217  inline EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_rhs.cols(); }
218 
219  template<typename Dest> inline void evalTo(Dest& dst) const
220  {
221  if(!is_same_dense(dst,m_rhs))
222  dst = m_rhs;
223  m_triangularMatrix.template solveInPlace<Side>(dst);
224  }
225 
226  protected:
227  const TriangularType& m_triangularMatrix;
228  typename Rhs::Nested m_rhs;
229 };
230 
231 } // namespace internal
232 
233 } // end namespace Eigen
234 
235 #endif // EIGEN_SOLVETRIANGULAR_H
SCALAR Scalar
Definition: bench_gemm.cpp:46
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
Expression of the transpose of a matrix.
Definition: Transpose.h:52
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
const unsigned int RowMajorBit
Definition: Constants.h:66
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Generic expression where a coefficient-wise binary operator is applied to two expressions.
Definition: CwiseBinaryOp.h:77
ReturnByValue< triangular_solve_retval > Base
triangular_solve_retval(const TriangularType &tri, const Rhs &rhs)
#define EIGEN_NOEXCEPT
Definition: Macros.h:1418
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
#define eigen_assert(x)
Definition: Macros.h:1037
#define EIGEN_CONSTEXPR
Definition: Macros.h:787
EIGEN_DEVICE_FUNC bool is_same_dense(const T1 &mat1, const T2 &mat2, typename enable_if< possibly_same_dense< T1, T2 >::value >::type *=0)
Definition: XprHelper.h:695
#define ei_declare_aligned_stack_constructed_variable(TYPE, NAME, SIZE, BUFFER)
Definition: Memory.h:768
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
EIGEN_CONSTEXPR Index size(const T &x)
Definition: Meta.h:479
#define EIGEN_DEVICE_FUNC
Definition: Macros.h:976
Tridiagonalization< MatrixXf > tri
remove_all< typename Rhs::Nested >::type RhsNestedCleaned
int EIGEN_BLAS_FUNC() copy(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
Definition: level1_impl.h:29
const int Dynamic
Definition: Constants.h:22
Generic expression where a coefficient-wise unary operator is applied to an expression.
Definition: CwiseUnaryOp.h:55
internal::plain_matrix_type_column_major< Rhs >::type ReturnType
The matrix class, also used for vectors and row-vectors.
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Definition: pytypes.h:1370


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:53