Home | History | Annotate | Download | only in IterativeSolvers 
      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 (at) 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
     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 
     94 namespace internal {
     95 
     96 template<typename _MatrixType, typename Rhs>
     97 struct 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