1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_INCOMPLETE_LU_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_INCOMPLETE_LU_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename _Scalar> 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass IncompleteLU 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _Scalar Scalar; 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> Vector; 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Vector::Index Index; 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef SparseMatrix<Scalar,RowMajor> FactorType; 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IncompleteLU() : m_isInitialized(false) {} 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IncompleteLU(const MatrixType& mat) : m_isInitialized(false) 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(mat); 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows() const { return m_lu.rows(); } 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols() const { return m_lu.cols(); } 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IncompleteLU& compute(const MatrixType& mat) 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_lu = mat; 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int size = mat.cols(); 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Vector diag(size); 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<size; ++i) 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename FactorType::InnerIterator k_it(m_lu,i); 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(; k_it && k_it.index()<i; ++k_it) 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int k = k_it.index(); 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath k_it.valueRef() /= diag(k); 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename FactorType::InnerIterator j_it(k_it); 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename FactorType::InnerIterator kj_it(m_lu, k); 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while(kj_it && kj_it.index()<=k) ++kj_it; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(++j_it; j_it; ) 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(kj_it.index()==j_it.index()) 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath j_it.valueRef() -= k_it.value() * kj_it.value(); 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++j_it; 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++kj_it; 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(kj_it.index()<j_it.index()) ++kj_it; 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else ++j_it; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(k_it && k_it.index()==i) diag(i) = k_it.value(); 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else diag(i) = 1; 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs, typename Dest> 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void _solve(const Rhs& b, Dest& x) const 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x = m_lu.template triangularView<UnitLower>().solve(b); 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x = m_lu.template triangularView<Upper>().solve(x); 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs> inline const internal::solve_retval<IncompleteLU, Rhs> 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solve(const MatrixBase<Rhs>& b) const 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "IncompleteLU is not initialized."); 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(cols()==b.rows() 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && "IncompleteLU::solve(): invalid number of rows of the right hand side matrix b"); 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return internal::solve_retval<IncompleteLU, Rhs>(*this, b.derived()); 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath FactorType m_lu; 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, typename Rhs> 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct solve_retval<IncompleteLU<_MatrixType>, Rhs> 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : solve_retval_base<IncompleteLU<_MatrixType>, Rhs> 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef IncompleteLU<_MatrixType> Dec; 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void evalTo(Dest& dst) const 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dec()._solve(rhs(),dst); 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_INCOMPLETE_LU_H 114