10 #ifndef EIGEN_SPARSETRIANGULARSOLVER_H 11 #define EIGEN_SPARSETRIANGULARSOLVER_H 17 template<
typename Lhs,
typename Rhs,
int Mode,
18 int UpLo = (Mode &
Lower)
23 int StorageOrder = int(traits<Lhs>::Flags) &
RowMajorBit>
24 struct sparse_solve_triangular_selector;
27 template<
typename Lhs,
typename Rhs,
int Mode>
28 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,
Lower,
RowMajor>
30 typedef typename Rhs::Scalar Scalar;
31 typedef evaluator<Lhs> LhsEval;
32 typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
33 static void run(
const Lhs& lhs, Rhs& other)
36 for(
Index col=0 ; col<other.cols() ; ++col)
38 for(
Index i=0; i<lhs.rows(); ++i)
40 Scalar tmp = other.coeff(i,col);
43 for(LhsIterator it(lhsEval, i); it; ++it)
46 lastIndex = it.index();
49 tmp -= lastVal * other.coeff(lastIndex,col);
52 other.coeffRef(i,col) = tmp;
55 eigen_assert(lastIndex==i);
56 other.coeffRef(i,col) = tmp/lastVal;
64 template<
typename Lhs,
typename Rhs,
int Mode>
65 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,
RowMajor>
67 typedef typename Rhs::Scalar Scalar;
68 typedef evaluator<Lhs> LhsEval;
69 typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
70 static void run(
const Lhs& lhs, Rhs& other)
73 for(
Index col=0 ; col<other.cols() ; ++col)
75 for(
Index i=lhs.rows()-1 ; i>=0 ; --i)
77 Scalar tmp = other.coeff(i,col);
79 LhsIterator it(lhsEval, i);
80 while(it && it.index()<i)
84 eigen_assert(it && it.index()==i);
88 else if (it && it.index() == i)
92 tmp -= it.value() * other.coeff(it.index(),col);
95 if (Mode & UnitDiag) other.coeffRef(i,col) = tmp;
96 else other.coeffRef(i,col) = tmp/l_ii;
103 template<
typename Lhs,
typename Rhs,
int Mode>
104 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,
Lower,
ColMajor>
106 typedef typename Rhs::Scalar Scalar;
107 typedef evaluator<Lhs> LhsEval;
108 typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
109 static void run(
const Lhs& lhs, Rhs& other)
111 LhsEval lhsEval(lhs);
112 for(
Index col=0 ; col<other.cols() ; ++col)
114 for(
Index i=0; i<lhs.cols(); ++i)
116 Scalar& tmp = other.coeffRef(i,col);
119 LhsIterator it(lhsEval, i);
120 while(it && it.index()<i)
122 if(!(Mode & UnitDiag))
124 eigen_assert(it && it.index()==i);
127 if (it && it.index()==i)
130 other.coeffRef(it.index(), col) -= tmp * it.value();
138 template<
typename Lhs,
typename Rhs,
int Mode>
139 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,
ColMajor>
141 typedef typename Rhs::Scalar Scalar;
142 typedef evaluator<Lhs> LhsEval;
143 typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
144 static void run(
const Lhs& lhs, Rhs& other)
146 LhsEval lhsEval(lhs);
147 for(
Index col=0 ; col<other.cols() ; ++col)
149 for(
Index i=lhs.cols()-1; i>=0; --i)
151 Scalar& tmp = other.coeffRef(i,col);
154 if(!(Mode & UnitDiag))
157 LhsIterator it(lhsEval, i);
158 while(it && it.index()!=i)
160 eigen_assert(it && it.index()==i);
161 other.coeffRef(i,col) /= it.value();
163 LhsIterator it(lhsEval, i);
164 for(; it && it.index()<i; ++it)
165 other.coeffRef(it.index(), col) -= tmp * it.value();
174 #ifndef EIGEN_PARSED_BY_DOXYGEN 176 template<
typename ExpressionType,
unsigned int Mode>
177 template<
typename OtherDerived>
178 void TriangularViewImpl<ExpressionType,Mode,Sparse>::solveInPlace(MatrixBase<OtherDerived>& other)
const 180 eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
181 eigen_assert((!(Mode &
ZeroDiag)) &&
bool(Mode & (Upper|
Lower)));
183 enum { copy = internal::traits<OtherDerived>::Flags &
RowMajorBit };
185 typedef typename internal::conditional<copy,
186 typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
187 OtherCopy otherCopy(other.derived());
189 internal::sparse_solve_triangular_selector<ExpressionType, typename internal::remove_reference<OtherCopy>::type, Mode>::run(derived().nestedExpression(), otherCopy);
200 template<
typename Lhs,
typename Rhs,
int Mode,
201 int UpLo = (Mode &
Lower)
206 int StorageOrder = int(Lhs::Flags) & (
RowMajorBit)>
207 struct sparse_solve_triangular_sparse_selector;
210 template<
typename Lhs,
typename Rhs,
int Mode,
int UpLo>
211 struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
213 typedef typename Rhs::Scalar Scalar;
214 typedef typename promote_index_type<typename traits<Lhs>::StorageIndex,
215 typename traits<Rhs>::StorageIndex>::type StorageIndex;
216 static void run(
const Lhs& lhs, Rhs& other)
218 const bool IsLower = (UpLo==
Lower);
219 AmbiVector<Scalar,StorageIndex> tempVector(other.rows()*2);
220 tempVector.setBounds(0,other.rows());
222 Rhs res(other.rows(), other.cols());
223 res.reserve(other.nonZeros());
225 for(
Index col=0 ; col<other.cols() ; ++col)
228 tempVector.init(.99);
229 tempVector.setZero();
230 tempVector.restart();
231 for (
typename Rhs::InnerIterator rhsIt(other, col); rhsIt; ++rhsIt)
233 tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
236 for(
Index i=IsLower?0:lhs.cols()-1;
237 IsLower?i<lhs.cols():i>=0;
240 tempVector.restart();
241 Scalar& ci = tempVector.coeffRef(i);
245 typename Lhs::InnerIterator it(lhs, i);
246 if(!(Mode & UnitDiag))
250 eigen_assert(it.index()==i);
254 ci /= lhs.coeff(i,i);
256 tempVector.restart();
262 tempVector.coeffRef(it.index()) -= ci * it.value();
266 for(; it && it.index()<i; ++it)
267 tempVector.coeffRef(it.index()) -= ci * it.value();
275 for (
typename AmbiVector<Scalar,StorageIndex>::Iterator it(tempVector); it; ++it)
281 res.insert(it.index(), col) = it.value();
286 other = res.markAsRValue();
292 #ifndef EIGEN_PARSED_BY_DOXYGEN 293 template<
typename ExpressionType,
unsigned int Mode>
294 template<
typename OtherDerived>
295 void TriangularViewImpl<ExpressionType,Mode,Sparse>::solveInPlace(SparseMatrixBase<OtherDerived>& other)
const 297 eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
298 eigen_assert( (!(Mode & ZeroDiag)) &&
bool(Mode & (Upper|
Lower)));
306 internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(derived().nestedExpression(), other.derived());
315 #endif // EIGEN_SPARSETRIANGULARSOLVER_H Definition: Constants.h:319
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:211
Definition: Constants.h:215
Definition: Eigen_Colamd.h:50
Definition: Constants.h:321