25 #ifndef EIGEN_SPARSETRIANGULARSOLVER_H
26 #define EIGEN_SPARSETRIANGULARSOLVER_H
32 template<
typename Lhs,
typename Rhs,
int Mode,
33 int UpLo = (Mode &
Lower)
39 struct sparse_solve_triangular_selector;
42 template<
typename Lhs,
typename Rhs,
int Mode>
43 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,
Lower,
RowMajor>
45 typedef typename Rhs::Scalar Scalar;
46 static void run(
const Lhs& lhs, Rhs& other)
50 for(
int i=0; i<lhs.rows(); ++i)
52 Scalar tmp = other.coeff(i,
col);
55 for(
typename Lhs::InnerIterator it(lhs, i); it; ++it)
58 lastIndex = it.index();
61 tmp -= lastVal * other.coeff(lastIndex,
col);
64 other.coeffRef(i,
col) = tmp;
68 other.coeffRef(i,
col) = tmp/lastVal;
76 template<
typename Lhs,
typename Rhs,
int Mode>
77 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,
RowMajor>
79 typedef typename Rhs::Scalar Scalar;
80 static void run(
const Lhs& lhs, Rhs& other)
84 for(
int i=lhs.rows()-1 ; i>=0 ; --i)
86 Scalar tmp = other.coeff(i,
col);
88 typename Lhs::InnerIterator it(lhs, i);
89 while(it && it.index()<i)
97 else if (it && it.index() == i)
101 tmp -= it.value() * other.coeff(it.index(),
col);
105 other.coeffRef(i,
col) = tmp;
107 other.coeffRef(i,
col) = tmp/l_ii;
114 template<
typename Lhs,
typename Rhs,
int Mode>
115 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,
Lower,
ColMajor>
117 typedef typename Rhs::Scalar Scalar;
118 static void run(
const Lhs& lhs, Rhs& other)
120 for(
int col=0 ;
col<other.cols() ; ++
col)
122 for(
int i=0; i<lhs.cols(); ++i)
124 Scalar& tmp = other.coeffRef(i,
col);
127 typename Lhs::InnerIterator it(lhs, i);
128 while(it && it.index()<i)
130 if(!(Mode & UnitDiag))
135 if (it && it.index()==i)
138 other.coeffRef(it.index(),
col) -= tmp * it.value();
146 template<
typename Lhs,
typename Rhs,
int Mode>
147 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,
ColMajor>
149 typedef typename Rhs::Scalar Scalar;
150 static void run(
const Lhs& lhs, Rhs& other)
152 for(
int col=0 ;
col<other.cols() ; ++
col)
154 for(
int i=lhs.cols()-1; i>=0; --i)
156 Scalar& tmp = other.coeffRef(i,
col);
159 if(!(Mode & UnitDiag))
162 typename Lhs::ReverseInnerIterator it(lhs, i);
163 while(it && it.index()!=i)
166 other.coeffRef(i,
col) /= it.value();
168 typename Lhs::InnerIterator it(lhs, i);
169 for(; it && it.index()<i; ++it)
170 other.coeffRef(it.index(),
col) -= tmp * it.value();
179 template<
typename ExpressionType,
int Mode>
180 template<
typename OtherDerived>
183 eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows());
186 enum { copy = internal::traits<OtherDerived>::Flags &
RowMajorBit };
188 typedef typename internal::conditional<copy,
190 OtherCopy otherCopy(other.derived());
192 internal::sparse_solve_triangular_selector<ExpressionType, typename internal::remove_reference<OtherCopy>::type, Mode>::run(m_matrix, otherCopy);
198 template<
typename ExpressionType,
int Mode>
199 template<
typename OtherDerived>
212 template<
typename Lhs,
typename Rhs,
int Mode,
213 int UpLo = (Mode &
Lower)
219 struct sparse_solve_triangular_sparse_selector;
222 template<
typename Lhs,
typename Rhs,
int Mode,
int UpLo>
223 struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
225 typedef typename Rhs::Scalar Scalar;
226 typedef typename promote_index_type<typename traits<Lhs>::Index,
228 static void run(
const Lhs& lhs, Rhs& other)
230 const bool IsLower = (UpLo==
Lower);
231 AmbiVector<Scalar,Index> tempVector(other.rows()*2);
232 tempVector.setBounds(0,other.rows());
234 Rhs res(other.rows(), other.cols());
235 res.reserve(other.nonZeros());
237 for(
int col=0 ;
col<other.cols() ; ++
col)
240 tempVector.init(.99);
241 tempVector.setZero();
242 tempVector.restart();
243 for (
typename Rhs::InnerIterator rhsIt(other,
col); rhsIt; ++rhsIt)
245 tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
248 for(
int i=IsLower?0:lhs.cols()-1;
249 IsLower?i<lhs.cols():i>=0;
252 tempVector.restart();
253 Scalar& ci = tempVector.coeffRef(i);
257 typename Lhs::InnerIterator it(lhs, i);
258 if(!(Mode & UnitDiag))
266 ci /= lhs.coeff(i,i);
268 tempVector.restart();
274 tempVector.coeffRef(it.index()) -= ci * it.value();
278 for(; it && it.index()<i; ++it)
279 tempVector.coeffRef(it.index()) -= ci * it.value();
287 for (
typename AmbiVector<Scalar,Index>::Iterator it(tempVector); it; ++it)
293 res.insert(it.index(),
col) = it.value();
298 other = res.markAsRValue();
304 template<
typename ExpressionType,
int Mode>
305 template<
typename OtherDerived>
308 eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.
rows());
317 internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(m_matrix, other.
derived());
323 #ifdef EIGEN2_SUPPORT
328 template<
typename Derived>
329 template<
typename OtherDerived>
332 this->
template triangular<Flags&(Upper|Lower)>().solveInPlace(other);
336 template<
typename Derived>
337 template<
typename OtherDerived>
338 typename internal::plain_matrix_type_column_major<OtherDerived>::type
339 SparseMatrixBase<Derived>::solveTriangular(
const MatrixBase<OtherDerived>& other)
const
341 typename internal::plain_matrix_type_column_major<OtherDerived>::type res(other);
342 derived().solveTriangularInPlace(res);
345 #endif // EIGEN2_SUPPORT
349 #endif // EIGEN_SPARSETRIANGULARSOLVER_H