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 : 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 30 IncompleteLU() {} 31 32 template<typename MatrixType> 33 IncompleteLU(const MatrixType& mat) 34 { 35 compute(mat); 36 } 37 38 Index rows() const { return m_lu.rows(); } 39 Index cols() const { return m_lu.cols(); } 40 41 template<typename MatrixType> 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> 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
