$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 00017 { 00018 typedef _Scalar Scalar; 00019 typedef Matrix<Scalar,Dynamic,1> Vector; 00020 typedef typename Vector::Index Index; 00021 typedef SparseMatrix<Scalar,RowMajor> FactorType; 00022 00023 public: 00024 typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; 00025 00026 IncompleteLU() : m_isInitialized(false) {} 00027 00028 template<typename MatrixType> 00029 IncompleteLU(const MatrixType& mat) : m_isInitialized(false) 00030 { 00031 compute(mat); 00032 } 00033 00034 Index rows() const { return m_lu.rows(); } 00035 Index cols() const { return m_lu.cols(); } 00036 00037 template<typename MatrixType> 00038 IncompleteLU& compute(const MatrixType& mat) 00039 { 00040 m_lu = mat; 00041 int size = mat.cols(); 00042 Vector diag(size); 00043 for(int i=0; i<size; ++i) 00044 { 00045 typename FactorType::InnerIterator k_it(m_lu,i); 00046 for(; k_it && k_it.index()<i; ++k_it) 00047 { 00048 int k = k_it.index(); 00049 k_it.valueRef() /= diag(k); 00050 00051 typename FactorType::InnerIterator j_it(k_it); 00052 typename FactorType::InnerIterator kj_it(m_lu, k); 00053 while(kj_it && kj_it.index()<=k) ++kj_it; 00054 for(++j_it; j_it; ) 00055 { 00056 if(kj_it.index()==j_it.index()) 00057 { 00058 j_it.valueRef() -= k_it.value() * kj_it.value(); 00059 ++j_it; 00060 ++kj_it; 00061 } 00062 else if(kj_it.index()<j_it.index()) ++kj_it; 00063 else ++j_it; 00064 } 00065 } 00066 if(k_it && k_it.index()==i) diag(i) = k_it.value(); 00067 else diag(i) = 1; 00068 } 00069 m_isInitialized = true; 00070 return *this; 00071 } 00072 00073 template<typename Rhs, typename Dest> 00074 void _solve(const Rhs& b, Dest& x) const 00075 { 00076 x = m_lu.template triangularView<UnitLower>().solve(b); 00077 x = m_lu.template triangularView<Upper>().solve(x); 00078 } 00079 00080 template<typename Rhs> inline const internal::solve_retval<IncompleteLU, Rhs> 00081 solve(const MatrixBase<Rhs>& b) const 00082 { 00083 eigen_assert(m_isInitialized && "IncompleteLU is not initialized."); 00084 eigen_assert(cols()==b.rows() 00085 && "IncompleteLU::solve(): invalid number of rows of the right hand side matrix b"); 00086 return internal::solve_retval<IncompleteLU, Rhs>(*this, b.derived()); 00087 } 00088 00089 protected: 00090 FactorType m_lu; 00091 bool m_isInitialized; 00092 }; 00093 00094 namespace internal { 00095 00096 template<typename _MatrixType, typename Rhs> 00097 struct solve_retval<IncompleteLU<_MatrixType>, Rhs> 00098 : solve_retval_base<IncompleteLU<_MatrixType>, Rhs> 00099 { 00100 typedef IncompleteLU<_MatrixType> Dec; 00101 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 00102 00103 template<typename Dest> void evalTo(Dest& dst) const 00104 { 00105 dec()._solve(rhs(),dst); 00106 } 00107 }; 00108 00109 } // end namespace internal 00110 00111 } // end namespace Eigen 00112 00113 #endif // EIGEN_INCOMPLETE_LU_H