IncompleteLU.h
Go to the documentation of this file.
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 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_INCOMPLETE_LU_H
11 #define EIGEN_INCOMPLETE_LU_H
12 
13 namespace Eigen {
14 
15 template <typename _Scalar>
17 {
18  typedef _Scalar Scalar;
20  typedef typename Vector::Index Index;
22 
23  public:
25 
27 
28  template<typename MatrixType>
29  IncompleteLU(const MatrixType& mat) : m_isInitialized(false)
30  {
31  compute(mat);
32  }
33 
34  Index rows() const { return m_lu.rows(); }
35  Index cols() const { return m_lu.cols(); }
36 
37  template<typename MatrixType>
38  IncompleteLU& compute(const MatrixType& mat)
39  {
40  m_lu = mat;
41  int size = mat.cols();
42  Vector diag(size);
43  for(int i=0; i<size; ++i)
44  {
45  typename FactorType::InnerIterator k_it(m_lu,i);
46  for(; k_it && k_it.index()<i; ++k_it)
47  {
48  int k = k_it.index();
49  k_it.valueRef() /= diag(k);
50 
51  typename FactorType::InnerIterator j_it(k_it);
52  typename FactorType::InnerIterator kj_it(m_lu, k);
53  while(kj_it && kj_it.index()<=k) ++kj_it;
54  for(++j_it; j_it; )
55  {
56  if(kj_it.index()==j_it.index())
57  {
58  j_it.valueRef() -= k_it.value() * kj_it.value();
59  ++j_it;
60  ++kj_it;
61  }
62  else if(kj_it.index()<j_it.index()) ++kj_it;
63  else ++j_it;
64  }
65  }
66  if(k_it && k_it.index()==i) diag(i) = k_it.value();
67  else diag(i) = 1;
68  }
69  m_isInitialized = true;
70  return *this;
71  }
72 
73  template<typename Rhs, typename Dest>
74  void _solve(const Rhs& b, Dest& x) const
75  {
76  x = m_lu.template triangularView<UnitLower>().solve(b);
77  x = m_lu.template triangularView<Upper>().solve(x);
78  }
79 
80  template<typename Rhs> inline const internal::solve_retval<IncompleteLU, Rhs>
81  solve(const MatrixBase<Rhs>& b) const
82  {
83  eigen_assert(m_isInitialized && "IncompleteLU is not initialized.");
84  eigen_assert(cols()==b.rows()
85  && "IncompleteLU::solve(): invalid number of rows of the right hand side matrix b");
86  return internal::solve_retval<IncompleteLU, Rhs>(*this, b.derived());
87  }
88 
89  protected:
90  FactorType m_lu;
92 };
93 
94 namespace internal {
95 
96 template<typename _MatrixType, typename Rhs>
97 struct solve_retval<IncompleteLU<_MatrixType>, Rhs>
98  : solve_retval_base<IncompleteLU<_MatrixType>, Rhs>
99 {
101  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
102 
103  template<typename Dest> void evalTo(Dest& dst) const
104  {
105  dec()._solve(rhs(),dst);
106  }
107 };
108 
109 } // end namespace internal
110 
111 } // end namespace Eigen
112 
113 #endif // EIGEN_INCOMPLETE_LU_H
Matrix< Scalar, Dynamic, 1 > Vector
Definition: IncompleteLU.h:19
Vector::Index Index
Definition: IncompleteLU.h:20
Index cols() const
Definition: SparseMatrix.h:121
Index cols() const
Definition: IncompleteLU.h:35
internal::traits< Derived >::Index Index
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
Matrix< Scalar, Dynamic, Dynamic > MatrixType
Definition: IncompleteLU.h:24
const internal::solve_retval< IncompleteLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: IncompleteLU.h:81
IncompleteLU & compute(const MatrixType &mat)
Definition: IncompleteLU.h:38
void rhs(const real_t *x, real_t *f)
void _solve(const Rhs &b, Dest &x) const
Definition: IncompleteLU.h:74
IncompleteLU(const MatrixType &mat)
Definition: IncompleteLU.h:29
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
Definition: Solve.h:61
SparseMatrix< Scalar, RowMajor > FactorType
Definition: IncompleteLU.h:21
EIGEN_STRONG_INLINE Index cols() const
#define eigen_assert(x)
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Index rows() const
Definition: SparseMatrix.h:119
Index rows() const
Definition: IncompleteLU.h:34


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:41