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 template<
typename ExpressionType,
unsigned int Mode>
175 template<
typename OtherDerived>
176 void TriangularViewImpl<ExpressionType,Mode,Sparse>::solveInPlace(MatrixBase<OtherDerived>& other)
const 178 eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
179 eigen_assert((!(Mode &
ZeroDiag)) &&
bool(Mode & (Upper|
Lower)));
181 enum { copy = internal::traits<OtherDerived>::Flags &
RowMajorBit };
183 typedef typename internal::conditional<copy,
184 typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
185 OtherCopy otherCopy(other.derived());
187 internal::sparse_solve_triangular_selector<ExpressionType, typename internal::remove_reference<OtherCopy>::type, Mode>::run(derived().nestedExpression(), otherCopy);
197 template<
typename Lhs,
typename Rhs,
int Mode,
198 int UpLo = (Mode &
Lower)
203 int StorageOrder = int(Lhs::Flags) & (
RowMajorBit)>
204 struct sparse_solve_triangular_sparse_selector;
207 template<
typename Lhs,
typename Rhs,
int Mode,
int UpLo>
208 struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
210 typedef typename Rhs::Scalar Scalar;
211 typedef typename promote_index_type<typename traits<Lhs>::StorageIndex,
212 typename traits<Rhs>::StorageIndex>::type StorageIndex;
213 static void run(
const Lhs& lhs, Rhs& other)
215 const bool IsLower = (UpLo==
Lower);
216 AmbiVector<Scalar,StorageIndex> tempVector(other.rows()*2);
217 tempVector.setBounds(0,other.rows());
219 Rhs res(other.rows(), other.cols());
220 res.reserve(other.nonZeros());
222 for(
Index col=0 ; col<other.cols() ; ++col)
225 tempVector.init(.99);
226 tempVector.setZero();
227 tempVector.restart();
228 for (
typename Rhs::InnerIterator rhsIt(other, col); rhsIt; ++rhsIt)
230 tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
233 for(
Index i=IsLower?0:lhs.cols()-1;
234 IsLower?i<lhs.cols():i>=0;
237 tempVector.restart();
238 Scalar& ci = tempVector.coeffRef(i);
242 typename Lhs::InnerIterator it(lhs, i);
243 if(!(Mode & UnitDiag))
247 eigen_assert(it.index()==i);
251 ci /= lhs.coeff(i,i);
253 tempVector.restart();
259 tempVector.coeffRef(it.index()) -= ci * it.value();
263 for(; it && it.index()<i; ++it)
264 tempVector.coeffRef(it.index()) -= ci * it.value();
272 for (
typename AmbiVector<Scalar,StorageIndex>::Iterator it(tempVector); it; ++it)
278 res.insert(it.index(), col) = it.value();
283 other = res.markAsRValue();
289 template<
typename ExpressionType,
unsigned int Mode>
290 template<
typename OtherDerived>
291 void TriangularViewImpl<ExpressionType,Mode,Sparse>::solveInPlace(SparseMatrixBase<OtherDerived>& other)
const 293 eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
294 eigen_assert( (!(Mode & ZeroDiag)) &&
bool(Mode & (Upper|
Lower)));
302 internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(derived().nestedExpression(), other.derived());
310 #endif // EIGEN_SPARSETRIANGULARSOLVER_H Definition: Constants.h:320
Namespace containing all symbols from the Eigen library.
Definition: Core:271
const unsigned int RowMajorBit
Definition: Constants.h:61
Definition: Constants.h:204
Definition: Constants.h:208
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: XprHelper.h:35
Definition: Constants.h:206
Definition: Constants.h:210
Definition: Eigen_Colamd.h:50
Definition: Constants.h:322