IncompleteLU.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_INCOMPLETE_LU_H
26 #define EIGEN_INCOMPLETE_LU_H
27 
28 namespace Eigen {
29 
30 template <typename _Scalar>
31 class IncompleteLU
32 {
33  typedef _Scalar Scalar;
34  typedef Matrix<Scalar,Dynamic,1> Vector;
35  typedef typename Vector::Index Index;
36  typedef SparseMatrix<Scalar,RowMajor> FactorType;
37 
38  public:
39  typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
40 
41  IncompleteLU() : m_isInitialized(false) {}
42 
43  template<typename MatrixType>
44  IncompleteLU(const MatrixType& mat) : m_isInitialized(false)
45  {
46  compute(mat);
47  }
48 
49  Index rows() const { return m_lu.rows(); }
50  Index cols() const { return m_lu.cols(); }
51 
52  template<typename MatrixType>
53  IncompleteLU& compute(const MatrixType& mat)
54  {
55  m_lu = mat;
56  int size = mat.cols();
57  Vector diag(size);
58  for(int i=0; i<size; ++i)
59  {
60  typename FactorType::InnerIterator k_it(m_lu,i);
61  for(; k_it && k_it.index()<i; ++k_it)
62  {
63  int k = k_it.index();
64  k_it.valueRef() /= diag(k);
65 
66  typename FactorType::InnerIterator j_it(k_it);
67  typename FactorType::InnerIterator kj_it(m_lu, k);
68  while(kj_it && kj_it.index()<=k) ++kj_it;
69  for(++j_it; j_it; )
70  {
71  if(kj_it.index()==j_it.index())
72  {
73  j_it.valueRef() -= k_it.value() * kj_it.value();
74  ++j_it;
75  ++kj_it;
76  }
77  else if(kj_it.index()<j_it.index()) ++kj_it;
78  else ++j_it;
79  }
80  }
81  if(k_it && k_it.index()==i) diag(i) = k_it.value();
82  else diag(i) = 1;
83  }
84  m_isInitialized = true;
85  return *this;
86  }
87 
88  template<typename Rhs, typename Dest>
89  void _solve(const Rhs& b, Dest& x) const
90  {
91  x = m_lu.template triangularView<UnitLower>().solve(b);
92  x = m_lu.template triangularView<Upper>().solve(x);
93  }
94 
95  template<typename Rhs> inline const internal::solve_retval<IncompleteLU, Rhs>
96  solve(const MatrixBase<Rhs>& b) const
97  {
98  eigen_assert(m_isInitialized && "IncompleteLU is not initialized.");
99  eigen_assert(cols()==b.rows()
100  && "IncompleteLU::solve(): invalid number of rows of the right hand side matrix b");
101  return internal::solve_retval<IncompleteLU, Rhs>(*this, b.derived());
102  }
103 
104  protected:
105  FactorType m_lu;
106  bool m_isInitialized;
107 };
108 
109 namespace internal {
110 
111 template<typename _MatrixType, typename Rhs>
112 struct solve_retval<IncompleteLU<_MatrixType>, Rhs>
113  : solve_retval_base<IncompleteLU<_MatrixType>, Rhs>
114 {
115  typedef IncompleteLU<_MatrixType> Dec;
116  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
117 
118  template<typename Dest> void evalTo(Dest& dst) const
119  {
120  dec()._solve(rhs(),dst);
121  }
122 };
123 
124 } // end namespace internal
125 
126 } // end namespace Eigen
127 
128 #endif // EIGEN_INCOMPLETE_LU_H