Eigen  3.3.4
SolveTriangular.h
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>
20 struct triangular_solve_vector;
21 
22 template <typename Scalar, typename Index, int Side, int Mode, bool Conjugate, int TriStorageOrder, int OtherStorageOrder>
23 struct triangular_solve_matrix;
24 
25 // small helper struct extracting some traits on the underlying solver operation
26 template<typename Lhs, typename Rhs, int Side>
27 class trsolve_traits
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)
36  ? CompleteUnrolling : NoUnrolling,
37  RhsVectors = RhsIsVectorAtCompileTime ? 1 : Dynamic
38  };
39 };
40 
41 template<typename Lhs, typename Rhs,
42  int Side, // can be OnTheLeft/OnTheRight
43  int Mode, // can be Upper/Lower | UnitDiag
44  int Unrolling = trsolve_traits<Lhs,Rhs,Side>::Unrolling,
45  int RhsVectors = trsolve_traits<Lhs,Rhs,Side>::RhsVectors
46  >
47 struct triangular_solver_selector;
48 
49 template<typename Lhs, typename Rhs, int Side, int Mode>
50 struct triangular_solver_selector<Lhs,Rhs,Side,Mode,NoUnrolling,1>
51 {
52  typedef typename Lhs::Scalar LhsScalar;
53  typedef typename Rhs::Scalar RhsScalar;
54  typedef blas_traits<Lhs> LhsProductTraits;
55  typedef typename LhsProductTraits::ExtractType ActualLhsType;
56  typedef Map<Matrix<RhsScalar,Dynamic,1>, Aligned> MappedRhs;
57  static 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>
82 struct triangular_solver_selector<Lhs,Rhs,Side,Mode,NoUnrolling,Dynamic>
83 {
84  typedef typename Rhs::Scalar Scalar;
85  typedef blas_traits<Lhs> LhsProductTraits;
86  typedef typename LhsProductTraits::DirectLinearAccessType ActualLhsType;
87 
88  static 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>
102  ::run(size, othersize, &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &rhs.coeffRef(0,0), 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>
112 struct triangular_solver_unroller;
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 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 
130  triangular_solver_unroller<Lhs,Rhs,Mode,LoopIndex+1,Size>::run(lhs,rhs);
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 void run(const Lhs&, Rhs&) {}
137 };
138 
139 template<typename Lhs, typename Rhs, int Mode>
140 struct triangular_solver_selector<Lhs,Rhs,OnTheLeft,Mode,CompleteUnrolling,1> {
141  static void run(const Lhs& lhs, Rhs& rhs)
142  { triangular_solver_unroller<Lhs,Rhs,Mode,0,Rhs::SizeAtCompileTime>::run(lhs,rhs); }
143 };
144 
145 template<typename Lhs, typename Rhs, int Mode>
146 struct triangular_solver_selector<Lhs,Rhs,OnTheRight,Mode,CompleteUnrolling,1> {
147  static void run(const Lhs& lhs, Rhs& rhs)
148  {
149  Transpose<const Lhs> trLhs(lhs);
150  Transpose<Rhs> trRhs(rhs);
151 
152  triangular_solver_unroller<Transpose<const Lhs>,Transpose<Rhs>,
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>
167 void TriangularViewImpl<MatrixType,Mode,Dense>::solveInPlace(const MatrixBase<OtherDerived>& _other) const
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((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
172 
173  enum { copy = (internal::traits<OtherDerived>::Flags & RowMajorBit) && OtherDerived::IsVectorAtCompileTime && OtherDerived::SizeAtCompileTime!=1};
174  typedef typename internal::conditional<copy,
175  typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
176  OtherCopy otherCopy(other);
177 
178  internal::triangular_solver_selector<MatrixType, typename internal::remove_reference<OtherCopy>::type,
179  Side, Mode>::run(derived().nestedExpression(), otherCopy);
180 
181  if (copy)
182  other = otherCopy;
183 }
184 
185 template<typename Derived, unsigned int Mode>
186 template<int Side, typename Other>
187 const internal::triangular_solve_retval<Side,TriangularView<Derived,Mode>,Other>
188 TriangularViewImpl<Derived,Mode,Dense>::solve(const MatrixBase<Other>& other) const
189 {
190  return internal::triangular_solve_retval<Side,TriangularViewType,Other>(derived(), other.derived());
191 }
192 #endif
193 
194 namespace internal {
195 
196 
197 template<int Side, typename TriangularType, typename Rhs>
198 struct traits<triangular_solve_retval<Side, TriangularType, Rhs> >
199 {
200  typedef typename internal::plain_matrix_type_column_major<Rhs>::type ReturnType;
201 };
202 
203 template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval
204  : public ReturnByValue<triangular_solve_retval<Side, TriangularType, Rhs> >
205 {
206  typedef typename remove_all<typename Rhs::Nested>::type RhsNestedCleaned;
207  typedef ReturnByValue<triangular_solve_retval> Base;
208 
209  triangular_solve_retval(const TriangularType& tri, const Rhs& rhs)
210  : m_triangularMatrix(tri), m_rhs(rhs)
211  {}
212 
213  inline Index rows() const { return m_rhs.rows(); }
214  inline Index cols() const { return m_rhs.cols(); }
215 
216  template<typename Dest> inline void evalTo(Dest& dst) const
217  {
218  if(!is_same_dense(dst,m_rhs))
219  dst = m_rhs;
220  m_triangularMatrix.template solveInPlace<Side>(dst);
221  }
222 
223  protected:
224  const TriangularType& m_triangularMatrix;
225  typename Rhs::Nested m_rhs;
226 };
227 
228 } // namespace internal
229 
230 } // end namespace Eigen
231 
232 #endif // EIGEN_SOLVETRIANGULAR_H
Definition: Constants.h:320
Definition: Constants.h:335
Namespace containing all symbols from the Eigen library.
Definition: Core:287
const unsigned int RowMajorBit
Definition: Constants.h:61
Definition: Constants.h:204
Definition: Constants.h:208
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Definition: Constants.h:235
Definition: Constants.h:333
Definition: Constants.h:206
Definition: Constants.h:210
Definition: Eigen_Colamd.h:50
Definition: Constants.h:322
const int Dynamic
Definition: Constants.h:21