GMRES.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 // Copyright (C) 2012, 2014 Kolja Brix <brix@igpm.rwth-aaachen.de>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_GMRES_H
12 #define EIGEN_GMRES_H
13 
14 namespace Eigen {
15 
16 namespace internal {
17 
55 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
56 bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Preconditioner & precond,
57  Index &iters, const Index &restart, typename Dest::RealScalar & tol_error) {
58 
59  using std::sqrt;
60  using std::abs;
61 
62  typedef typename Dest::RealScalar RealScalar;
63  typedef typename Dest::Scalar Scalar;
66 
67  const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
68 
69  if(rhs.norm() <= considerAsZero)
70  {
71  x.setZero();
72  tol_error = 0;
73  return true;
74  }
75 
76  RealScalar tol = tol_error;
77  const Index maxIters = iters;
78  iters = 0;
79 
80  const Index m = mat.rows();
81 
82  // residual and preconditioned residual
83  VectorType p0 = rhs - mat*x;
84  VectorType r0 = precond.solve(p0);
85 
86  const RealScalar r0Norm = r0.norm();
87 
88  // is initial guess already good enough?
89  if(r0Norm == 0)
90  {
91  tol_error = 0;
92  return true;
93  }
94 
95  // storage for Hessenberg matrix and Householder data
96  FMatrixType H = FMatrixType::Zero(m, restart + 1);
97  VectorType w = VectorType::Zero(restart + 1);
98  VectorType tau = VectorType::Zero(restart + 1);
99 
100  // storage for Jacobi rotations
101  std::vector < JacobiRotation < Scalar > > G(restart);
102 
103  // storage for temporaries
104  VectorType t(m), v(m), workspace(m), x_new(m);
105 
106  // generate first Householder vector
107  Ref<VectorType> H0_tail = H.col(0).tail(m - 1);
108  RealScalar beta;
109  r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta);
110  w(0) = Scalar(beta);
111 
112  for (Index k = 1; k <= restart; ++k)
113  {
114  ++iters;
115 
116  v = VectorType::Unit(m, k - 1);
117 
118  // apply Householder reflections H_{1} ... H_{k-1} to v
119  // TODO: use a HouseholderSequence
120  for (Index i = k - 1; i >= 0; --i) {
121  v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
122  }
123 
124  // apply matrix M to v: v = mat * v;
125  t.noalias() = mat * v;
126  v = precond.solve(t);
127 
128  // apply Householder reflections H_{k-1} ... H_{1} to v
129  // TODO: use a HouseholderSequence
130  for (Index i = 0; i < k; ++i) {
131  v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
132  }
133 
134  if (v.tail(m - k).norm() != 0.0)
135  {
136  if (k <= restart)
137  {
138  // generate new Householder vector
139  Ref<VectorType> Hk_tail = H.col(k).tail(m - k - 1);
140  v.tail(m - k).makeHouseholder(Hk_tail, tau.coeffRef(k), beta);
141 
142  // apply Householder reflection H_{k} to v
143  v.tail(m - k).applyHouseholderOnTheLeft(Hk_tail, tau.coeffRef(k), workspace.data());
144  }
145  }
146 
147  if (k > 1)
148  {
149  for (Index i = 0; i < k - 1; ++i)
150  {
151  // apply old Givens rotations to v
152  v.applyOnTheLeft(i, i + 1, G[i].adjoint());
153  }
154  }
155 
156  if (k<m && v(k) != (Scalar) 0)
157  {
158  // determine next Givens rotation
159  G[k - 1].makeGivens(v(k - 1), v(k));
160 
161  // apply Givens rotation to v and w
162  v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
163  w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
164  }
165 
166  // insert coefficients into upper matrix triangle
167  H.col(k-1).head(k) = v.head(k);
168 
169  tol_error = abs(w(k)) / r0Norm;
170  bool stop = (k==m || tol_error < tol || iters == maxIters);
171 
172  if (stop || k == restart)
173  {
174  // solve upper triangular system
175  Ref<VectorType> y = w.head(k);
176  H.topLeftCorner(k, k).template triangularView <Upper>().solveInPlace(y);
177 
178  // use Horner-like scheme to calculate solution vector
179  x_new.setZero();
180  for (Index i = k - 1; i >= 0; --i)
181  {
182  x_new(i) += y(i);
183  // apply Householder reflection H_{i} to x_new
184  x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
185  }
186 
187  x += x_new;
188 
189  if(stop)
190  {
191  return true;
192  }
193  else
194  {
195  k=0;
196 
197  // reset data for restart
198  p0.noalias() = rhs - mat*x;
199  r0 = precond.solve(p0);
200 
201  // clear Hessenberg matrix and Householder data
202  H.setZero();
203  w.setZero();
204  tau.setZero();
205 
206  // generate first Householder vector
207  r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta);
208  w(0) = Scalar(beta);
209  }
210  }
211  }
212 
213  return false;
214 
215 }
216 
217 }
218 
219 template< typename _MatrixType,
220  typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
221 class GMRES;
222 
223 namespace internal {
224 
225 template< typename _MatrixType, typename _Preconditioner>
226 struct traits<GMRES<_MatrixType,_Preconditioner> >
227 {
228  typedef _MatrixType MatrixType;
229  typedef _Preconditioner Preconditioner;
230 };
231 
232 }
233 
268 template< typename _MatrixType, typename _Preconditioner>
269 class GMRES : public IterativeSolverBase<GMRES<_MatrixType,_Preconditioner> >
270 {
272  using Base::matrix;
273  using Base::m_error;
274  using Base::m_iterations;
275  using Base::m_info;
276  using Base::m_isInitialized;
277 
278 private:
280 
281 public:
282  using Base::_solve_impl;
283  typedef _MatrixType MatrixType;
284  typedef typename MatrixType::Scalar Scalar;
286  typedef _Preconditioner Preconditioner;
287 
288 public:
289 
291  GMRES() : Base(), m_restart(30) {}
292 
303  template<typename MatrixDerived>
304  explicit GMRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_restart(30) {}
305 
306  ~GMRES() {}
307 
310  Index get_restart() { return m_restart; }
311 
315  void set_restart(const Index restart) { m_restart=restart; }
316 
318  template<typename Rhs,typename Dest>
319  void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
320  {
321  m_iterations = Base::maxIterations();
322  m_error = Base::m_tolerance;
323  bool ret = internal::gmres(matrix(), b, x, Base::m_preconditioner, m_iterations, m_restart, m_error);
324  m_info = (!ret) ? NumericalIssue
325  : m_error <= Base::m_tolerance ? Success
326  : NoConvergence;
327  }
328 
329 protected:
330 
331 };
332 
333 } // end namespace Eigen
334 
335 #endif // EIGEN_GMRES_H
Index m_restart
Definition: GMRES.h:279
Matrix3f m
A preconditioner based on the digonal entries.
SCALAR Scalar
Definition: bench_gemm.cpp:46
MatrixType::RealScalar RealScalar
Definition: GMRES.h:285
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
Scalar * b
Definition: benchVecAdd.cpp:17
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
Index get_restart()
Definition: GMRES.h:310
JacobiRotation< float > G
#define min(a, b)
Definition: datatypes.h:19
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
MatrixXf MatrixType
Vector3f p0
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics set mxtics default set mytics default set mx2tics default set my2tics default set xtics border mirror norotate autofreq set ytics border mirror norotate autofreq set ztics border nomirror norotate autofreq set nox2tics set noy2tics set timestamp bottom norotate set rrange [*:*] noreverse nowriteback set trange [*:*] noreverse nowriteback set urange [*:*] noreverse nowriteback set vrange [*:*] noreverse nowriteback set xlabel matrix size set x2label set timefmt d m y n H
_MatrixType MatrixType
Definition: GMRES.h:283
void _solve_vector_with_guess_impl(const Rhs &b, Dest &x) const
Definition: GMRES.h:319
_Preconditioner Preconditioner
Definition: GMRES.h:286
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Array< int, Dynamic, 1 > v
GMRES(const EigenBase< MatrixDerived > &A)
Definition: GMRES.h:304
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
RowVector3d w
DenseIndex ret
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:281
void set_restart(const Index restart)
Definition: GMRES.h:315
IterativeSolverBase< GMRES > Base
Definition: GMRES.h:271
A GMRES solver for sparse square problems.
Definition: GMRES.h:221
MatrixType::Scalar Scalar
Definition: GMRES.h:284
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
const G double tol
Definition: Group.h:86
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
#define abs(x)
Definition: datatypes.h:17
Base class for linear iterative solvers.
Point2 t(10, 10)
bool gmres(const MatrixType &mat, const Rhs &rhs, Dest &x, const Preconditioner &precond, Index &iters, const Index &restart, typename Dest::RealScalar &tol_error)
Definition: GMRES.h:56


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:19