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