Eigen  3.2.91
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 Index, int Size,
111  bool Stop = Index==Size>
112 struct triangular_solver_unroller;
113 
114 template<typename Lhs, typename Rhs, int Mode, int Index, int Size>
115 struct triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,false> {
116  enum {
117  IsLower = ((Mode&Lower)==Lower),
118  I = IsLower ? Index : Size - Index - 1,
119  S = IsLower ? 0 : I+1
120  };
121  static void run(const Lhs& lhs, Rhs& rhs)
122  {
123  if (Index>0)
124  rhs.coeffRef(I) -= lhs.row(I).template segment<Index>(S).transpose()
125  .cwiseProduct(rhs.template segment<Index>(S)).sum();
126 
127  if(!(Mode & UnitDiag))
128  rhs.coeffRef(I) /= lhs.coeff(I,I);
129 
130  triangular_solver_unroller<Lhs,Rhs,Mode,Index+1,Size>::run(lhs,rhs);
131  }
132 };
133 
134 template<typename Lhs, typename Rhs, int Mode, int Index, int Size>
135 struct triangular_solver_unroller<Lhs,Rhs,Mode,Index,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 
171 template<typename MatrixType, unsigned int Mode>
172 template<int Side, typename OtherDerived>
173 void TriangularViewImpl<MatrixType,Mode,Dense>::solveInPlace(const MatrixBase<OtherDerived>& _other) const
174 {
175  OtherDerived& other = _other.const_cast_derived();
176  eigen_assert( derived().cols() == derived().rows() && ((Side==OnTheLeft && derived().cols() == other.rows()) || (Side==OnTheRight && derived().cols() == other.cols())) );
177  eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
178 
179  enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit && OtherDerived::IsVectorAtCompileTime };
180  typedef typename internal::conditional<copy,
181  typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
182  OtherCopy otherCopy(other);
183 
184  internal::triangular_solver_selector<MatrixType, typename internal::remove_reference<OtherCopy>::type,
185  Side, Mode>::run(derived().nestedExpression(), otherCopy);
186 
187  if (copy)
188  other = otherCopy;
189 }
190 
212 template<typename Derived, unsigned int Mode>
213 template<int Side, typename Other>
214 const internal::triangular_solve_retval<Side,TriangularView<Derived,Mode>,Other>
215 TriangularViewImpl<Derived,Mode,Dense>::solve(const MatrixBase<Other>& other) const
216 {
217  return internal::triangular_solve_retval<Side,TriangularViewType,Other>(derived(), other.derived());
218 }
219 
220 namespace internal {
221 
222 
223 template<int Side, typename TriangularType, typename Rhs>
224 struct traits<triangular_solve_retval<Side, TriangularType, Rhs> >
225 {
226  typedef typename internal::plain_matrix_type_column_major<Rhs>::type ReturnType;
227 };
228 
229 template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval
230  : public ReturnByValue<triangular_solve_retval<Side, TriangularType, Rhs> >
231 {
232  typedef typename remove_all<typename Rhs::Nested>::type RhsNestedCleaned;
233  typedef ReturnByValue<triangular_solve_retval> Base;
234 
235  triangular_solve_retval(const TriangularType& tri, const Rhs& rhs)
236  : m_triangularMatrix(tri), m_rhs(rhs)
237  {}
238 
239  inline Index rows() const { return m_rhs.rows(); }
240  inline Index cols() const { return m_rhs.cols(); }
241 
242  template<typename Dest> inline void evalTo(Dest& dst) const
243  {
244  if(!(is_same<RhsNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_rhs)))
245  dst = m_rhs;
246  m_triangularMatrix.template solveInPlace<Side>(dst);
247  }
248 
249  protected:
250  const TriangularType& m_triangularMatrix;
251  typename Rhs::Nested m_rhs;
252 };
253 
254 } // namespace internal
255 
256 } // end namespace Eigen
257 
258 #endif // EIGEN_SOLVETRIANGULAR_H
Definition: Constants.h:314
Definition: Constants.h:196
Definition: Constants.h:202
Definition: LDLT.h:16
Definition: Constants.h:327
const unsigned int RowMajorBit
Definition: Constants.h:53
Definition: Constants.h:198
Definition: Constants.h:200
Definition: Eigen_Colamd.h:54
Definition: Constants.h:312
Definition: Constants.h:325
Definition: Constants.h:227