| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| #ifndef EIGEN_INCOMPLETE_LU_H |
| #define EIGEN_INCOMPLETE_LU_H |
| |
| namespace Eigen { |
| |
| template <typename _Scalar> |
| class IncompleteLU : public SparseSolverBase<IncompleteLU<_Scalar> > |
| { |
| protected: |
| typedef SparseSolverBase<IncompleteLU<_Scalar> > Base; |
| using Base::m_isInitialized; |
| |
| typedef _Scalar Scalar; |
| typedef Matrix<Scalar,Dynamic,1> Vector; |
| typedef typename Vector::Index Index; |
| typedef SparseMatrix<Scalar,RowMajor> FactorType; |
| |
| public: |
| typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; |
| |
| IncompleteLU() {} |
| |
| template<typename MatrixType> |
| IncompleteLU(const MatrixType& mat) |
| { |
| compute(mat); |
| } |
| |
| Index rows() const { return m_lu.rows(); } |
| Index cols() const { return m_lu.cols(); } |
| |
| template<typename MatrixType> |
| IncompleteLU& compute(const MatrixType& mat) |
| { |
| m_lu = mat; |
| int size = mat.cols(); |
| Vector diag(size); |
| for(int i=0; i<size; ++i) |
| { |
| typename FactorType::InnerIterator k_it(m_lu,i); |
| for(; k_it && k_it.index()<i; ++k_it) |
| { |
| int k = k_it.index(); |
| k_it.valueRef() /= diag(k); |
| |
| typename FactorType::InnerIterator j_it(k_it); |
| typename FactorType::InnerIterator kj_it(m_lu, k); |
| while(kj_it && kj_it.index()<=k) ++kj_it; |
| for(++j_it; j_it; ) |
| { |
| if(kj_it.index()==j_it.index()) |
| { |
| j_it.valueRef() -= k_it.value() * kj_it.value(); |
| ++j_it; |
| ++kj_it; |
| } |
| else if(kj_it.index()<j_it.index()) ++kj_it; |
| else ++j_it; |
| } |
| } |
| if(k_it && k_it.index()==i) diag(i) = k_it.value(); |
| else diag(i) = 1; |
| } |
| m_isInitialized = true; |
| return *this; |
| } |
| |
| template<typename Rhs, typename Dest> |
| void _solve_impl(const Rhs& b, Dest& x) const |
| { |
| x = m_lu.template triangularView<UnitLower>().solve(b); |
| x = m_lu.template triangularView<Upper>().solve(x); |
| } |
| |
| protected: |
| FactorType m_lu; |
| }; |
| |
| } // end namespace Eigen |
| |
| #endif // EIGEN_INCOMPLETE_LU_H |