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
17{
18    typedef _Scalar Scalar;
19    typedef Matrix<Scalar,Dynamic,1> Vector;
20    typedef typename Vector::Index Index;
21    typedef SparseMatrix<Scalar,RowMajor> FactorType;
22
23  public:
24    typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
25
26    IncompleteLU() : m_isInitialized(false) {}
27
28    template<typename MatrixType>
29    IncompleteLU(const MatrixType& mat) : m_isInitialized(false)
30    {
31      compute(mat);
32    }
33
34    Index rows() const { return m_lu.rows(); }
35    Index cols() const { return m_lu.cols(); }
36
37    template<typename MatrixType>
38    IncompleteLU& compute(const MatrixType& mat)
39    {
40      m_lu = mat;
41      int size = mat.cols();
42      Vector diag(size);
43      for(int i=0; i<size; ++i)
44      {
45        typename FactorType::InnerIterator k_it(m_lu,i);
46        for(; k_it && k_it.index()<i; ++k_it)
47        {
48          int k = k_it.index();
49          k_it.valueRef() /= diag(k);
50
51          typename FactorType::InnerIterator j_it(k_it);
52          typename FactorType::InnerIterator kj_it(m_lu, k);
53          while(kj_it && kj_it.index()<=k) ++kj_it;
54          for(++j_it; j_it; )
55          {
56            if(kj_it.index()==j_it.index())
57            {
58              j_it.valueRef() -= k_it.value() * kj_it.value();
59              ++j_it;
60              ++kj_it;
61            }
62            else if(kj_it.index()<j_it.index()) ++kj_it;
63            else                                ++j_it;
64          }
65        }
66        if(k_it && k_it.index()==i) diag(i) = k_it.value();
67        else                        diag(i) = 1;
68      }
69      m_isInitialized = true;
70      return *this;
71    }
72
73    template<typename Rhs, typename Dest>
74    void _solve(const Rhs& b, Dest& x) const
75    {
76      x = m_lu.template triangularView<UnitLower>().solve(b);
77      x = m_lu.template triangularView<Upper>().solve(x);
78    }
79
80    template<typename Rhs> inline const internal::solve_retval<IncompleteLU, Rhs>
81    solve(const MatrixBase<Rhs>& b) const
82    {
83      eigen_assert(m_isInitialized && "IncompleteLU is not initialized.");
84      eigen_assert(cols()==b.rows()
85                && "IncompleteLU::solve(): invalid number of rows of the right hand side matrix b");
86      return internal::solve_retval<IncompleteLU, Rhs>(*this, b.derived());
87    }
88
89  protected:
90    FactorType m_lu;
91    bool m_isInitialized;
92};
93
94namespace internal {
95
96template<typename _MatrixType, typename Rhs>
97struct solve_retval<IncompleteLU<_MatrixType>, Rhs>
98  : solve_retval_base<IncompleteLU<_MatrixType>, Rhs>
99{
100  typedef IncompleteLU<_MatrixType> Dec;
101  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
102
103  template<typename Dest> void evalTo(Dest& dst) const
104  {
105    dec()._solve(rhs(),dst);
106  }
107};
108
109} // end namespace internal
110
111} // end namespace Eigen
112
113#endif // EIGEN_INCOMPLETE_LU_H
114