10 #ifndef EIGEN_SOLVETRIANGULAR_H 11 #define EIGEN_SOLVETRIANGULAR_H 19 template<
typename LhsScalar,
typename RhsScalar,
typename Index,
int S
ide,
int Mode,
bool Conjugate,
int StorageOrder>
20 struct triangular_solve_vector;
22 template <
typename Scalar,
typename Index,
int S
ide,
int Mode,
bool Conjugate,
int TriStorageOrder,
int OtherStorageOrder,
int OtherInnerStr
ide>
23 struct triangular_solve_matrix;
26 template<
typename Lhs,
typename Rhs,
int S
ide>
31 RhsIsVectorAtCompileTime = (Side==
OnTheLeft ? Rhs::ColsAtCompileTime : Rhs::RowsAtCompileTime)==1
35 Unrolling = (RhsIsVectorAtCompileTime && Rhs::SizeAtCompileTime !=
Dynamic && Rhs::SizeAtCompileTime <= 8)
36 ? CompleteUnrolling : NoUnrolling,
37 RhsVectors = RhsIsVectorAtCompileTime ? 1 :
Dynamic 41 template<
typename Lhs,
typename Rhs,
44 int Unrolling = trsolve_traits<Lhs,Rhs,Side>::Unrolling,
45 int RhsVectors = trsolve_traits<Lhs,Rhs,Side>::RhsVectors
47 struct triangular_solver_selector;
49 template<
typename Lhs,
typename Rhs,
int S
ide,
int Mode>
50 struct triangular_solver_selector<Lhs,Rhs,Side,Mode,NoUnrolling,1>
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 EIGEN_DEVICE_FUNC
void run(
const Lhs& lhs, Rhs& rhs)
59 ActualLhsType actualLhs = LhsProductTraits::extract(lhs);
63 bool useRhsDirectly = Rhs::InnerStrideAtCompileTime==1 || rhs.innerStride()==1;
65 ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhs,rhs.size(),
66 (useRhsDirectly ? rhs.data() : 0));
69 MappedRhs(actualRhs,rhs.size()) = rhs;
71 triangular_solve_vector<LhsScalar, RhsScalar,
Index, Side, Mode, LhsProductTraits::NeedToConjugate,
73 ::run(actualLhs.cols(), actualLhs.data(), actualLhs.outerStride(), actualRhs);
76 rhs = MappedRhs(actualRhs, rhs.size());
81 template<
typename Lhs,
typename Rhs,
int S
ide,
int Mode>
82 struct triangular_solver_selector<Lhs,Rhs,Side,Mode,NoUnrolling,
Dynamic>
84 typedef typename Rhs::Scalar Scalar;
85 typedef blas_traits<Lhs> LhsProductTraits;
86 typedef typename LhsProductTraits::DirectLinearAccessType ActualLhsType;
88 static EIGEN_DEVICE_FUNC
void run(
const Lhs& lhs, Rhs& rhs)
90 typename internal::add_const_on_value_type<ActualLhsType>::type actualLhs = LhsProductTraits::extract(lhs);
92 const Index size = lhs.rows();
93 const Index othersize = Side==
OnTheLeft? rhs.cols() : rhs.rows();
96 Rhs::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime, Lhs::MaxRowsAtCompileTime,4> BlockingType;
98 BlockingType blocking(rhs.rows(), rhs.cols(), size, 1,
false);
100 triangular_solve_matrix<Scalar,
Index,Side,Mode,LhsProductTraits::NeedToConjugate,(int(Lhs::Flags) &
RowMajorBit) ?
RowMajor : ColMajor,
102 ::run(size, othersize, &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &rhs.coeffRef(0,0), rhs.innerStride(), rhs.outerStride(), blocking);
110 template<
typename Lhs,
typename Rhs,
int Mode,
int LoopIndex,
int Size,
111 bool Stop = LoopIndex==Size>
112 struct triangular_solver_unroller;
114 template<
typename Lhs,
typename Rhs,
int Mode,
int LoopIndex,
int Size>
115 struct triangular_solver_unroller<Lhs,Rhs,Mode,LoopIndex,Size,false> {
118 DiagIndex = IsLower ? LoopIndex : Size - LoopIndex - 1,
119 StartIndex = IsLower ? 0 : DiagIndex+1
121 static EIGEN_DEVICE_FUNC
void run(
const Lhs& lhs, Rhs& rhs)
124 rhs.coeffRef(DiagIndex) -= lhs.row(DiagIndex).template segment<LoopIndex>(StartIndex).transpose()
125 .cwiseProduct(rhs.template segment<LoopIndex>(StartIndex)).sum();
128 rhs.coeffRef(DiagIndex) /= lhs.coeff(DiagIndex,DiagIndex);
130 triangular_solver_unroller<Lhs,Rhs,Mode,LoopIndex+1,Size>::run(lhs,rhs);
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&) {}
139 template<
typename Lhs,
typename Rhs,
int Mode>
140 struct triangular_solver_selector<Lhs,Rhs,
OnTheLeft,Mode,CompleteUnrolling,1> {
141 static EIGEN_DEVICE_FUNC
void run(
const Lhs& lhs, Rhs& rhs)
142 { triangular_solver_unroller<Lhs,Rhs,Mode,0,Rhs::SizeAtCompileTime>::run(lhs,rhs); }
145 template<
typename Lhs,
typename Rhs,
int Mode>
146 struct triangular_solver_selector<Lhs,Rhs,
OnTheRight,Mode,CompleteUnrolling,1> {
147 static EIGEN_DEVICE_FUNC
void run(
const Lhs& lhs, Rhs& rhs)
149 Transpose<const Lhs> trLhs(lhs);
150 Transpose<Rhs> trRhs(rhs);
152 triangular_solver_unroller<Transpose<const Lhs>,Transpose<Rhs>,
154 0,Rhs::SizeAtCompileTime>::run(trLhs,trRhs);
164 #ifndef EIGEN_PARSED_BY_DOXYGEN 165 template<
typename MatrixType,
unsigned int Mode>
166 template<
int S
ide,
typename OtherDerived>
167 EIGEN_DEVICE_FUNC
void TriangularViewImpl<MatrixType,Mode,Dense>::solveInPlace(
const MatrixBase<OtherDerived>& _other)
const 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))));
173 if (derived().cols() == 0)
176 enum { copy = (internal::traits<OtherDerived>::Flags &
RowMajorBit) && OtherDerived::IsVectorAtCompileTime && OtherDerived::SizeAtCompileTime!=1};
177 typedef typename internal::conditional<copy,
178 typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
179 OtherCopy otherCopy(other);
181 internal::triangular_solver_selector<MatrixType, typename internal::remove_reference<OtherCopy>::type,
182 Side, Mode>::run(derived().nestedExpression(), otherCopy);
188 template<
typename Derived,
unsigned int Mode>
189 template<
int S
ide,
typename Other>
190 const internal::triangular_solve_retval<Side,TriangularView<Derived,Mode>,Other>
191 TriangularViewImpl<Derived,Mode,Dense>::solve(
const MatrixBase<Other>& other)
const 193 return internal::triangular_solve_retval<Side,TriangularViewType,Other>(derived(), other.derived());
200 template<
int S
ide,
typename TriangularType,
typename Rhs>
201 struct traits<triangular_solve_retval<Side, TriangularType, Rhs> >
203 typedef typename internal::plain_matrix_type_column_major<Rhs>::type ReturnType;
206 template<
int S
ide,
typename TriangularType,
typename Rhs>
struct triangular_solve_retval
207 :
public ReturnByValue<triangular_solve_retval<Side, TriangularType, Rhs> >
209 typedef typename remove_all<typename Rhs::Nested>::type RhsNestedCleaned;
210 typedef ReturnByValue<triangular_solve_retval> Base;
212 triangular_solve_retval(
const TriangularType& tri,
const Rhs& rhs)
213 : m_triangularMatrix(tri), m_rhs(rhs)
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(); }
219 template<
typename Dest>
inline void evalTo(Dest& dst)
const 221 if(!is_same_dense(dst,m_rhs))
223 m_triangularMatrix.template solveInPlace<Side>(dst);
227 const TriangularType& m_triangularMatrix;
228 typename Rhs::Nested m_rhs;
235 #endif // EIGEN_SOLVETRIANGULAR_H Definition: Constants.h:319
Definition: Constants.h:334
Namespace containing all symbols from the Eigen library.
Definition: Core:141
const unsigned int RowMajorBit
Definition: Constants.h:66
Definition: Constants.h:209
Definition: Constants.h:213
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: Constants.h:240
Definition: Constants.h:332
Definition: Constants.h:211
Definition: Constants.h:215
Definition: Eigen_Colamd.h:50
Definition: Constants.h:321
const int Dynamic
Definition: Constants.h:22