![]() |
Eigen-unsupported
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_INCOMPLETE_LU_H 00011 #define EIGEN_INCOMPLETE_LU_H 00012 00013 namespace Eigen { 00014 00015 template <typename _Scalar> 00016 class IncompleteLU : public SparseSolverBase<IncompleteLU<_Scalar> > 00017 { 00018 protected: 00019 typedef SparseSolverBase<IncompleteLU<_Scalar> > Base; 00020 using Base::m_isInitialized; 00021 00022 typedef _Scalar Scalar; 00023 typedef Matrix<Scalar,Dynamic,1> Vector; 00024 typedef typename Vector::Index Index; 00025 typedef SparseMatrix<Scalar,RowMajor> FactorType; 00026 00027 public: 00028 typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; 00029 00030 IncompleteLU() {} 00031 00032 template<typename MatrixType> 00033 IncompleteLU(const MatrixType& mat) 00034 { 00035 compute(mat); 00036 } 00037 00038 Index rows() const { return m_lu.rows(); } 00039 Index cols() const { return m_lu.cols(); } 00040 00041 template<typename MatrixType> 00042 IncompleteLU& compute(const MatrixType& mat) 00043 { 00044 m_lu = mat; 00045 int size = mat.cols(); 00046 Vector diag(size); 00047 for(int i=0; i<size; ++i) 00048 { 00049 typename FactorType::InnerIterator k_it(m_lu,i); 00050 for(; k_it && k_it.index()<i; ++k_it) 00051 { 00052 int k = k_it.index(); 00053 k_it.valueRef() /= diag(k); 00054 00055 typename FactorType::InnerIterator j_it(k_it); 00056 typename FactorType::InnerIterator kj_it(m_lu, k); 00057 while(kj_it && kj_it.index()<=k) ++kj_it; 00058 for(++j_it; j_it; ) 00059 { 00060 if(kj_it.index()==j_it.index()) 00061 { 00062 j_it.valueRef() -= k_it.value() * kj_it.value(); 00063 ++j_it; 00064 ++kj_it; 00065 } 00066 else if(kj_it.index()<j_it.index()) ++kj_it; 00067 else ++j_it; 00068 } 00069 } 00070 if(k_it && k_it.index()==i) diag(i) = k_it.value(); 00071 else diag(i) = 1; 00072 } 00073 m_isInitialized = true; 00074 return *this; 00075 } 00076 00077 template<typename Rhs, typename Dest> 00078 void _solve_impl(const Rhs& b, Dest& x) const 00079 { 00080 x = m_lu.template triangularView<UnitLower>().solve(b); 00081 x = m_lu.template triangularView<Upper>().solve(x); 00082 } 00083 00084 protected: 00085 FactorType m_lu; 00086 }; 00087 00088 } // end namespace Eigen 00089 00090 #endif // EIGEN_INCOMPLETE_LU_H