IncompleteLU.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_INCOMPLETE_LU_H
00011 #define EIGEN_INCOMPLETE_LU_H
00012 
00013 namespace Eigen { 
00014 
00015 template <typename _Scalar>
00016 class IncompleteLU : public SparseSolverBase<IncompleteLU<_Scalar> >
00017 {
00018   protected:
00019     typedef SparseSolverBase<IncompleteLU<_Scalar> > Base;
00020     using Base::m_isInitialized;
00021     
00022     typedef _Scalar Scalar;
00023     typedef Matrix<Scalar,Dynamic,1> Vector;
00024     typedef typename Vector::Index Index;
00025     typedef SparseMatrix<Scalar,RowMajor> FactorType;
00026 
00027   public:
00028     typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
00029 
00030     IncompleteLU() {}
00031 
00032     template<typename MatrixType>
00033     IncompleteLU(const MatrixType& mat)
00034     {
00035       compute(mat);
00036     }
00037 
00038     Index rows() const { return m_lu.rows(); }
00039     Index cols() const { return m_lu.cols(); }
00040 
00041     template<typename MatrixType>
00042     IncompleteLU& compute(const MatrixType& mat)
00043     {
00044       m_lu = mat;
00045       int size = mat.cols();
00046       Vector diag(size);
00047       for(int i=0; i<size; ++i)
00048       {
00049         typename FactorType::InnerIterator k_it(m_lu,i);
00050         for(; k_it && k_it.index()<i; ++k_it)
00051         {
00052           int k = k_it.index();
00053           k_it.valueRef() /= diag(k);
00054 
00055           typename FactorType::InnerIterator j_it(k_it);
00056           typename FactorType::InnerIterator kj_it(m_lu, k);
00057           while(kj_it && kj_it.index()<=k) ++kj_it;
00058           for(++j_it; j_it; )
00059           {
00060             if(kj_it.index()==j_it.index())
00061             {
00062               j_it.valueRef() -= k_it.value() * kj_it.value();
00063               ++j_it;
00064               ++kj_it;
00065             }
00066             else if(kj_it.index()<j_it.index()) ++kj_it;
00067             else                                ++j_it;
00068           }
00069         }
00070         if(k_it && k_it.index()==i) diag(i) = k_it.value();
00071         else                        diag(i) = 1;
00072       }
00073       m_isInitialized = true;
00074       return *this;
00075     }
00076 
00077     template<typename Rhs, typename Dest>
00078     void _solve_impl(const Rhs& b, Dest& x) const
00079     {
00080       x = m_lu.template triangularView<UnitLower>().solve(b);
00081       x = m_lu.template triangularView<Upper>().solve(x);
00082     }
00083 
00084   protected:
00085     FactorType m_lu;
00086 };
00087 
00088 } // end namespace Eigen
00089 
00090 #endif // EIGEN_INCOMPLETE_LU_H
 All Classes Functions Variables Typedefs Enumerator