25 #ifndef EIGEN_SOLVETRIANGULAR_H
26 #define EIGEN_SOLVETRIANGULAR_H
34 template<
typename LhsScalar,
typename RhsScalar,
typename Index,
int S
ide,
int Mode,
bool Conjugate,
int StorageOrder>
35 struct triangular_solve_vector;
37 template <
typename Scalar,
typename Index,
int S
ide,
int Mode,
bool Conjugate,
int TriStorageOrder,
int OtherStorageOrder>
38 struct triangular_solve_matrix;
41 template<
typename Lhs,
typename Rhs,
int S
ide>
46 RhsIsVectorAtCompileTime = (Side==
OnTheLeft ? Rhs::ColsAtCompileTime : Rhs::RowsAtCompileTime)==1
50 Unrolling = (RhsIsVectorAtCompileTime && Rhs::SizeAtCompileTime !=
Dynamic && Rhs::SizeAtCompileTime <= 8)
52 RhsVectors = RhsIsVectorAtCompileTime ? 1 :
Dynamic
56 template<
typename Lhs,
typename Rhs,
59 int Unrolling = trsolve_traits<Lhs,Rhs,Side>::Unrolling,
60 int RhsVectors = trsolve_traits<Lhs,Rhs,Side>::RhsVectors
62 struct triangular_solver_selector;
64 template<
typename Lhs,
typename Rhs,
int S
ide,
int Mode>
65 struct triangular_solver_selector<Lhs,Rhs,Side,Mode,
NoUnrolling,1>
67 typedef typename Lhs::Scalar LhsScalar;
68 typedef typename Rhs::Scalar RhsScalar;
69 typedef blas_traits<Lhs> LhsProductTraits;
70 typedef typename LhsProductTraits::ExtractType ActualLhsType;
71 typedef Map<Matrix<RhsScalar,Dynamic,1>,
Aligned> MappedRhs;
72 static void run(
const Lhs& lhs, Rhs& rhs)
74 ActualLhsType actualLhs = LhsProductTraits::extract(lhs);
78 bool useRhsDirectly = Rhs::InnerStrideAtCompileTime==1 || rhs.innerStride()==1;
81 (useRhsDirectly ? rhs.data() : 0));
84 MappedRhs(actualRhs,rhs.size()) = rhs;
86 triangular_solve_vector<LhsScalar, RhsScalar,
typename Lhs::Index, Side, Mode, LhsProductTraits::NeedToConjugate,
88 ::run(actualLhs.cols(), actualLhs.data(), actualLhs.outerStride(), actualRhs);
91 rhs = MappedRhs(actualRhs, rhs.size());
96 template<
typename Lhs,
typename Rhs,
int S
ide,
int Mode>
99 typedef typename Rhs::Scalar Scalar;
100 typedef typename Rhs::Index Index;
101 typedef blas_traits<Lhs> LhsProductTraits;
102 typedef typename LhsProductTraits::DirectLinearAccessType ActualLhsType;
103 static void run(
const Lhs& lhs, Rhs& rhs)
105 typename internal::add_const_on_value_type<ActualLhsType>::type actualLhs = LhsProductTraits::extract(lhs);
108 ::run(lhs.rows(), Side==
OnTheLeft? rhs.cols() : rhs.rows(), &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &rhs.coeffRef(0,0), rhs.outerStride());
116 template<
typename Lhs,
typename Rhs,
int Mode,
int Index,
int Size,
117 bool Stop = Index==Size>
118 struct triangular_solver_unroller;
120 template<
typename Lhs,
typename Rhs,
int Mode,
int Index,
int Size>
121 struct triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,false> {
124 I = IsLower ? Index : Size - Index - 1,
125 S = IsLower ? 0 : I+1
127 static void run(
const Lhs& lhs, Rhs& rhs)
130 rhs.coeffRef(I) -= lhs.row(I).template segment<Index>(S).transpose()
131 .cwiseProduct(rhs.template segment<Index>(S)).sum();
134 rhs.coeffRef(I) /= lhs.coeff(I,I);
136 triangular_solver_unroller<Lhs,Rhs,Mode,Index+1,Size>::run(lhs,rhs);
140 template<
typename Lhs,
typename Rhs,
int Mode,
int Index,
int Size>
141 struct triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,true> {
142 static void run(
const Lhs&, Rhs&) {}
145 template<
typename Lhs,
typename Rhs,
int Mode>
147 static void run(
const Lhs& lhs, Rhs& rhs)
148 { triangular_solver_unroller<Lhs,Rhs,Mode,0,Rhs::SizeAtCompileTime>::run(lhs,rhs); }
151 template<
typename Lhs,
typename Rhs,
int Mode>
153 static void run(
const Lhs& lhs, Rhs& rhs)
155 Transpose<const Lhs> trLhs(lhs);
156 Transpose<Rhs> trRhs(rhs);
158 triangular_solver_unroller<Transpose<const Lhs>,Transpose<Rhs>,
160 0,Rhs::SizeAtCompileTime>::run(trLhs,trRhs);
177 template<
typename MatrixType,
unsigned int Mode>
178 template<
int S
ide,
typename OtherDerived>
181 OtherDerived& other = _other.const_cast_derived();
185 enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit && OtherDerived::IsVectorAtCompileTime };
186 typedef typename internal::conditional<copy,
188 OtherCopy otherCopy(other);
190 internal::triangular_solver_selector<MatrixType, typename internal::remove_reference<OtherCopy>::type,
191 Side, Mode>::run(nestedExpression(), otherCopy);
218 template<
typename Derived,
unsigned int Mode>
219 template<
int S
ide,
typename Other>
220 const internal::triangular_solve_retval<Side,TriangularView<Derived,Mode>,Other>
223 return internal::triangular_solve_retval<Side,TriangularView,Other>(*
this, other.derived());
229 template<
int S
ide,
typename TriangularType,
typename Rhs>
230 struct traits<triangular_solve_retval<Side, TriangularType, Rhs> >
235 template<
int S
ide,
typename TriangularType,
typename Rhs>
struct triangular_solve_retval
236 :
public ReturnByValue<triangular_solve_retval<Side, TriangularType, Rhs> >
238 typedef typename remove_all<typename Rhs::Nested>::type RhsNestedCleaned;
239 typedef ReturnByValue<triangular_solve_retval> Base;
240 typedef typename Base::Index Index;
242 triangular_solve_retval(
const TriangularType& tri,
const Rhs& rhs)
243 : m_triangularMatrix(tri), m_rhs(rhs)
246 inline Index rows()
const {
return m_rhs.rows(); }
247 inline Index cols()
const {
return m_rhs.cols(); }
249 template<
typename Dest>
inline void evalTo(Dest& dst)
const
253 m_triangularMatrix.template solveInPlace<Side>(dst);
257 const TriangularType& m_triangularMatrix;
258 typename Rhs::Nested m_rhs;
265 #endif // EIGEN_SOLVETRIANGULAR_H