src/IterativeSolvers/IncompleteLU.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_INCOMPLETE_LU_H
0011 #define EIGEN_INCOMPLETE_LU_H
0012
0013 namespace Eigen {
0014
0015 template <typename _Scalar>
0016 class IncompleteLU : public SparseSolverBase<IncompleteLU<_Scalar> >
0017 {
0018   protected:
0019     typedef SparseSolverBase<IncompleteLU<_Scalar> > Base;
0020     using Base::m_isInitialized;
0021
0022     typedef _Scalar Scalar;
0023     typedef Matrix<Scalar,Dynamic,1> Vector;
0024     typedef typename Vector::Index Index;
0025     typedef SparseMatrix<Scalar,RowMajor> FactorType;
0026
0027   public:
0028     typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
0029
0030     IncompleteLU() {}
0031
0032     template<typename MatrixType>
0033     IncompleteLU(const MatrixType& mat)
0034     {
0035       compute(mat);
0036     }
0037
0038     Index rows() const { return m_lu.rows(); }
0039     Index cols() const { return m_lu.cols(); }
0040
0041     template<typename MatrixType>
0042     IncompleteLU& compute(const MatrixType& mat)
0043     {
0044       m_lu = mat;
0045       int size = mat.cols();
0046       Vector diag(size);
0047       for(int i=0; i<size; ++i)
0048       {
0049         typename FactorType::InnerIterator k_it(m_lu,i);
0050         for(; k_it && k_it.index()<i; ++k_it)
0051         {
0052           int k = k_it.index();
0053           k_it.valueRef() /= diag(k);
0054
0055           typename FactorType::InnerIterator j_it(k_it);
0056           typename FactorType::InnerIterator kj_it(m_lu, k);
0057           while(kj_it && kj_it.index()<=k) ++kj_it;
0058           for(++j_it; j_it; )
0059           {
0060             if(kj_it.index()==j_it.index())
0061             {
0062               j_it.valueRef() -= k_it.value() * kj_it.value();
0063               ++j_it;
0064               ++kj_it;
0065             }
0066             else if(kj_it.index()<j_it.index()) ++kj_it;
0067             else                                ++j_it;
0068           }
0069         }
0070         if(k_it && k_it.index()==i) diag(i) = k_it.value();
0071         else                        diag(i) = 1;
0072       }
0073       m_isInitialized = true;
0074       return *this;
0075     }
0076
0077     template<typename Rhs, typename Dest>
0078     void _solve_impl(const Rhs& b, Dest& x) const
0079     {
0080       x = m_lu.template triangularView<UnitLower>().solve(b);
0081       x = m_lu.template triangularView<Upper>().solve(x);
0082     }
0083
0084   protected:
0085     FactorType m_lu;
0086 };
0087
0088 } // end namespace Eigen
0089
0090 #endif // EIGEN_INCOMPLETE_LU_H