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 
13 namespace Eigen {
14 
15 template <typename _Scalar>
16 class 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 
IncompleteLU()30     IncompleteLU() {}
31 
32     template<typename MatrixType>
IncompleteLU(const MatrixType & mat)33     IncompleteLU(const MatrixType& mat)
34     {
35       compute(mat);
36     }
37 
rows()38     Index rows() const { return m_lu.rows(); }
cols()39     Index cols() const { return m_lu.cols(); }
40 
41     template<typename MatrixType>
compute(const MatrixType & mat)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>
_solve_impl(const Rhs & b,Dest & x)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