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 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 IncompleteLU()26 IncompleteLU() : m_isInitialized(false) {} 27 28 template<typename MatrixType> IncompleteLU(const MatrixType & mat)29 IncompleteLU(const MatrixType& mat) : m_isInitialized(false) 30 { 31 compute(mat); 32 } 33 rows()34 Index rows() const { return m_lu.rows(); } cols()35 Index cols() const { return m_lu.cols(); } 36 37 template<typename MatrixType> compute(const MatrixType & mat)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> _solve(const Rhs & b,Dest & x)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> solve(const MatrixBase<Rhs> & b)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