1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_INCOMPLETE_LU_H
11#define EIGEN_INCOMPLETE_LU_H
12
13namespace Eigen {
14
15template <typename _Scalar>
16class IncompleteLU : public SparseSolverBase<IncompleteLU<_Scalar> >
17{
18  protected:
19    typedef SparseSolverBase<IncompleteLU<_Scalar> > Base;
20    using Base::m_isInitialized;
21
22    typedef _Scalar Scalar;
23    typedef Matrix<Scalar,Dynamic,1> Vector;
24    typedef typename Vector::Index Index;
25    typedef SparseMatrix<Scalar,RowMajor> FactorType;
26
27  public:
28    typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
29
30    IncompleteLU() {}
31
32    template<typename MatrixType>
33    IncompleteLU(const MatrixType& mat)
34    {
35      compute(mat);
36    }
37
38    Index rows() const { return m_lu.rows(); }
39    Index cols() const { return m_lu.cols(); }
40
41    template<typename MatrixType>
42    IncompleteLU& compute(const MatrixType& mat)
43    {
44      m_lu = mat;
45      int size = mat.cols();
46      Vector diag(size);
47      for(int i=0; i<size; ++i)
48      {
49        typename FactorType::InnerIterator k_it(m_lu,i);
50        for(; k_it && k_it.index()<i; ++k_it)
51        {
52          int k = k_it.index();
53          k_it.valueRef() /= diag(k);
54
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;
58          for(++j_it; j_it; )
59          {
60            if(kj_it.index()==j_it.index())
61            {
62              j_it.valueRef() -= k_it.value() * kj_it.value();
63              ++j_it;
64              ++kj_it;
65            }
66            else if(kj_it.index()<j_it.index()) ++kj_it;
67            else                                ++j_it;
68          }
69        }
70        if(k_it && k_it.index()==i) diag(i) = k_it.value();
71        else                        diag(i) = 1;
72      }
73      m_isInitialized = true;
74      return *this;
75    }
76
77    template<typename Rhs, typename Dest>
78    void _solve_impl(const Rhs& b, Dest& x) const
79    {
80      x = m_lu.template triangularView<UnitLower>().solve(b);
81      x = m_lu.template triangularView<Upper>().solve(x);
82    }
83
84  protected:
85    FactorType m_lu;
86};
87
88} // end namespace Eigen
89
90#endif // EIGEN_INCOMPLETE_LU_H
91