10#ifndef EIGEN_INCOMPLETE_LU_H
11#define EIGEN_INCOMPLETE_LU_H
15template <
typename _Scalar>
16class IncompleteLU :
public SparseSolverBase<IncompleteLU<_Scalar> >
19 typedef SparseSolverBase<IncompleteLU<_Scalar> > Base;
20 using Base::m_isInitialized;
22 typedef _Scalar Scalar;
23 typedef Matrix<Scalar,Dynamic,1>
Vector;
24 typedef typename Vector::Index
Index;
25 typedef SparseMatrix<Scalar,RowMajor> FactorType;
28 typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
32 template<
typename MatrixType>
33 IncompleteLU(
const MatrixType& mat)
38 Index rows()
const {
return m_lu.rows(); }
39 Index cols()
const {
return m_lu.cols(); }
41 template<
typename MatrixType>
42 IncompleteLU& compute(
const MatrixType& mat)
45 int size = mat.cols();
47 for(
int i=0; i<size; ++i)
49 typename FactorType::InnerIterator k_it(m_lu,i);
50 for(; k_it && k_it.index()<i; ++k_it)
53 k_it.valueRef() /= diag(k);
55 typename FactorType::InnerIterator j_it(k_it);
56 typename FactorType::InnerIterator kj_it(m_lu, k);
57 while(kj_it && kj_it.index()<=k) ++kj_it;
60 if(kj_it.index()==j_it.index())
62 j_it.valueRef() -= k_it.value() * kj_it.value();
66 else if(kj_it.index()<j_it.index()) ++kj_it;
70 if(k_it && k_it.index()==i) diag(i) = k_it.value();
73 m_isInitialized =
true;
77 template<
typename Rhs,
typename Dest>
78 void _solve_impl(
const Rhs& b, Dest& x)
const
80 x = m_lu.template triangularView<UnitLower>().solve(b);
81 x = m_lu.template triangularView<Upper>().solve(x);
Matrix< Type, Size, 1 > Vector
Namespace containing all symbols from the Eigen library.
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index