25 #ifndef EIGEN_INCOMPLETE_LU_H
26 #define EIGEN_INCOMPLETE_LU_H
30 template <
typename _Scalar>
33 typedef _Scalar Scalar;
34 typedef Matrix<Scalar,Dynamic,1> Vector;
35 typedef typename Vector::Index Index;
36 typedef SparseMatrix<Scalar,RowMajor> FactorType;
39 typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
41 IncompleteLU() : m_isInitialized(false) {}
43 template<
typename MatrixType>
44 IncompleteLU(
const MatrixType& mat) : m_isInitialized(false)
49 Index rows()
const {
return m_lu.rows(); }
50 Index cols()
const {
return m_lu.cols(); }
52 template<
typename MatrixType>
53 IncompleteLU& compute(
const MatrixType& mat)
56 int size = mat.cols();
58 for(
int i=0; i<size; ++i)
60 typename FactorType::InnerIterator k_it(m_lu,i);
61 for(; k_it && k_it.index()<i; ++k_it)
64 k_it.valueRef() /= diag(k);
66 typename FactorType::InnerIterator j_it(k_it);
67 typename FactorType::InnerIterator kj_it(m_lu, k);
68 while(kj_it && kj_it.index()<=k) ++kj_it;
71 if(kj_it.index()==j_it.index())
73 j_it.valueRef() -= k_it.value() * kj_it.value();
77 else if(kj_it.index()<j_it.index()) ++kj_it;
81 if(k_it && k_it.index()==i) diag(i) = k_it.value();
84 m_isInitialized =
true;
88 template<
typename Rhs,
typename Dest>
89 void _solve(
const Rhs& b, Dest& x)
const
91 x = m_lu.template triangularView<UnitLower>().solve(b);
92 x = m_lu.template triangularView<Upper>().solve(x);
95 template<
typename Rhs>
inline const internal::solve_retval<IncompleteLU, Rhs>
96 solve(
const MatrixBase<Rhs>& b)
const
98 eigen_assert(m_isInitialized &&
"IncompleteLU is not initialized.");
99 eigen_assert(cols()==b.rows()
100 &&
"IncompleteLU::solve(): invalid number of rows of the right hand side matrix b");
101 return internal::solve_retval<IncompleteLU, Rhs>(*
this, b.derived());
106 bool m_isInitialized;
111 template<
typename _MatrixType,
typename Rhs>
112 struct solve_retval<IncompleteLU<_MatrixType>, Rhs>
113 : solve_retval_base<IncompleteLU<_MatrixType>, Rhs>
115 typedef IncompleteLU<_MatrixType> Dec;
116 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
118 template<typename Dest>
void evalTo(Dest& dst)
const
120 dec()._solve(rhs(),dst);
128 #endif // EIGEN_INCOMPLETE_LU_H