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;
104 static void run(
const Lhs& lhs, Rhs& rhs)
106 typename internal::add_const_on_value_type<ActualLhsType>::type actualLhs = LhsProductTraits::extract(lhs);
108 const Index size = lhs.rows();
109 const Index othersize = Side==
OnTheLeft? rhs.cols() : rhs.rows();
112 Rhs::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime, Lhs::MaxRowsAtCompileTime,4> BlockingType;
114 BlockingType blocking(rhs.rows(), rhs.cols(), size);
116 triangular_solve_matrix<Scalar,Index,Side,Mode,LhsProductTraits::NeedToConjugate,(
int(Lhs::Flags) &
RowMajorBit) ?
RowMajor : ColMajor,
118 ::run(size, othersize, &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &rhs.coeffRef(0,0), rhs.outerStride(), blocking);
126 template<
typename Lhs,
typename Rhs,
int Mode,
int Index,
int Size,
127 bool Stop = Index==Size>
128 struct triangular_solver_unroller;
130 template<
typename Lhs,
typename Rhs,
int Mode,
int Index,
int Size>
131 struct triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,false> {
134 I = IsLower ? Index : Size - Index - 1,
135 S = IsLower ? 0 : I+1
137 static void run(
const Lhs& lhs, Rhs& rhs)
140 rhs.coeffRef(I) -= lhs.row(I).template segment<Index>(S).transpose()
141 .cwiseProduct(rhs.template segment<Index>(S)).sum();
144 rhs.coeffRef(I) /= lhs.coeff(I,I);
146 triangular_solver_unroller<Lhs,Rhs,Mode,Index+1,Size>::run(lhs,rhs);
150 template<
typename Lhs,
typename Rhs,
int Mode,
int Index,
int Size>
151 struct triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,true> {
152 static void run(
const Lhs&, Rhs&) {}
155 template<
typename Lhs,
typename Rhs,
int Mode>
157 static void run(
const Lhs& lhs, Rhs& rhs)
158 { triangular_solver_unroller<Lhs,Rhs,Mode,0,Rhs::SizeAtCompileTime>::run(lhs,rhs); }
161 template<
typename Lhs,
typename Rhs,
int Mode>
163 static void run(
const Lhs& lhs, Rhs& rhs)
165 Transpose<const Lhs> trLhs(lhs);
166 Transpose<Rhs> trRhs(rhs);
168 triangular_solver_unroller<Transpose<const Lhs>,Transpose<Rhs>,
170 0,Rhs::SizeAtCompileTime>::run(trLhs,trRhs);
187 template<
typename MatrixType,
unsigned int Mode>
188 template<
int S
ide,
typename OtherDerived>
191 OtherDerived& other = _other.const_cast_derived();
195 enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit && OtherDerived::IsVectorAtCompileTime };
196 typedef typename internal::conditional<copy,
198 OtherCopy otherCopy(other);
200 internal::triangular_solver_selector<MatrixType, typename internal::remove_reference<OtherCopy>::type,
201 Side, Mode>::run(nestedExpression(), otherCopy);
228 template<
typename Derived,
unsigned int Mode>
229 template<
int S
ide,
typename Other>
230 const internal::triangular_solve_retval<Side,TriangularView<Derived,Mode>,Other>
233 return internal::triangular_solve_retval<Side,TriangularView,Other>(*
this, other.derived());
239 template<
int S
ide,
typename TriangularType,
typename Rhs>
240 struct traits<triangular_solve_retval<Side, TriangularType, Rhs> >
245 template<
int S
ide,
typename TriangularType,
typename Rhs>
struct triangular_solve_retval
246 :
public ReturnByValue<triangular_solve_retval<Side, TriangularType, Rhs> >
248 typedef typename remove_all<typename Rhs::Nested>::type RhsNestedCleaned;
249 typedef ReturnByValue<triangular_solve_retval> Base;
250 typedef typename Base::Index Index;
252 triangular_solve_retval(
const TriangularType& tri,
const Rhs& rhs)
253 : m_triangularMatrix(tri), m_rhs(rhs)
256 inline Index rows()
const {
return m_rhs.rows(); }
257 inline Index cols()
const {
return m_rhs.cols(); }
259 template<
typename Dest>
inline void evalTo(Dest& dst)
const
263 m_triangularMatrix.template solveInPlace<Side>(dst);
267 const TriangularType& m_triangularMatrix;
268 typename Rhs::Nested m_rhs;
275 #endif // EIGEN_SOLVETRIANGULAR_H