Block-Structured AMR Software Framework
AMReX_GMRES.H
Go to the documentation of this file.
1 #ifndef AMREX_GMRES_H_
2 #define AMREX_GMRES_H_
3 #include <AMReX_Config.H>
4 
5 #include <AMReX_BLProfiler.H>
6 #include <AMReX_Print.H>
7 #include <AMReX_TableData.H>
8 #include <AMReX_Vector.H>
9 #include <cmath>
10 #include <limits>
11 #include <memory>
12 
13 namespace amrex {
14 
59 template <typename V, typename M>
60 class GMRES
61 {
62 public:
63 
64  using RT = typename M::RT; // double or float
65 
66  GMRES ();
67 
71  void define (M& linop);
72 
82  void solve (V& a_sol, V const& a_rhs, RT a_tol_rel, RT a_tol_abs, int a_its=-1);
83 
85  void setVerbose (int v) { m_verbose = v; }
86 
88  void setRestartLength (int rl);
89 
91  void setMaxIters (int niters) { m_maxiter = niters; }
92 
94  [[nodiscard]] int getNumIters () const { return m_its; }
95 
97  [[nodiscard]] int getStatus () const { return m_status; }
98 
100  [[nodiscard]] RT getResidualNorm () const { return m_res; }
101 
102 private:
103  void clear ();
105  void cycle (V& a_xx, int& a_status, int& a_itcount, RT& a_rnorm0);
106  void build_solution (V& a_xx, int it);
107  void compute_residual (V& a_rr, V const& a_xx, V const& a_bb);
108 
109  bool converged (RT r0, RT r) const;
110 
112  void update_hessenberg (int it, bool happyend, RT& res);
113 
114  int m_verbose = 0;
115  int m_maxiter = 2000;
116  int m_its = 0;
117  int m_status = -1;
118  int m_restrtlen = 30;
120  RT m_rtol = RT(0);
121  RT m_atol = RT(0);
129  std::unique_ptr<V> m_v_tmp_rhs;
130  std::unique_ptr<V> m_v_tmp_lhs;
132  M* m_linop = nullptr;
133 };
134 
135 template <typename V, typename M>
137 {
138  allocate_scratch();
139 }
140 
141 template <typename V, typename M>
143 {
144  int rs = m_restrtlen;
145 
146  m_hh_1d.resize(std::size_t(rs + 2) * (rs + 1));
147  m_hh = Table2D<RT>(m_hh_1d.data(), {0,0}, {rs+1,rs}); // (0:rs+1,0:rs)
148 
149  m_hes_1d.resize(std::size_t(rs + 2) * (rs + 1));
150  m_hes = Table2D<RT>(m_hes_1d.data(), {0,0}, {rs+1,rs}); // (0:rs+1,0:rs)
151 
152  m_grs.resize(rs + 2);
153  m_cc.resize(rs + 1);
154  m_ss.resize(rs + 1);
155 }
156 
157 template <typename V, typename M>
159 {
160  if (m_restrtlen != rl) {
161  m_restrtlen = rl;
162  allocate_scratch();
163  m_vv.clear();
164  }
165 }
166 
167 template <typename V, typename M>
168 void GMRES<V,M>::define (M& linop)
169 {
170  clear();
171  m_linop = &linop;
172 }
173 
174 template <typename V, typename M>
176 {
177  m_its = 0;
178  m_status = -1;
180  m_v_tmp_rhs.reset();
181  m_v_tmp_lhs.reset();
182  m_vv.clear();
183  m_linop = nullptr;
184 }
185 
186 template <typename V, typename M>
187 bool GMRES<V,M>::converged (RT r0, RT r) const
188 {
189  return (r < r0*m_rtol) || (r < m_atol);
190 }
191 
192 template <typename V, typename M>
193 void GMRES<V,M>::solve (V& a_sol, V const& a_rhs, RT a_tol_rel, RT a_tol_abs, int a_its)
194 {
195  BL_PROFILE("GMRES::solve()");
196 
197  AMREX_ALWAYS_ASSERT(m_linop != nullptr);
198 
199  auto t0 = amrex::second();
200 
201  if (m_v_tmp_rhs == nullptr) {
202  m_v_tmp_rhs = std::make_unique<V>(m_linop->makeVecRHS());
203  }
204  if (m_v_tmp_lhs == nullptr) {
205  m_v_tmp_lhs = std::make_unique<V>(m_linop->makeVecLHS());
206  }
207  if (m_vv.empty()) {
208  m_vv.reserve(m_restrtlen+1);
209  for (int i = 0; i < 2; ++i) { // to save space, start with just 2
210  m_vv.emplace_back(m_linop->makeVecRHS());
211  }
212  }
213 
214  m_rtol = a_tol_rel;
215  m_atol = a_tol_abs;
216 
217  if (a_its < 0) { a_its = m_maxiter; }
218 
219  auto rnorm0 = RT(0);
220 
221  m_linop->assign(m_vv[0], a_rhs);
222  m_linop->setToZero(a_sol);
223 
224  m_its = 0;
225  m_status = -1;
226  cycle(a_sol, m_status, m_its, rnorm0);
227 
228  while (m_status == -1 && m_its < a_its) {
229  compute_residual(m_vv[0], a_sol, a_rhs);
230  cycle(a_sol, m_status, m_its, rnorm0);
231  }
232 
233  if (m_status == -1 && m_its >= a_its) { m_status = 1; }
234 
235  m_v_tmp_rhs.reset();
236  m_v_tmp_lhs.reset();
237  m_vv.clear();
238 
239  auto t1 = amrex::second();
240  if (m_verbose > 0) {
241  amrex::Print() << "GMRES: Solve Time = " << t1-t0 << '\n';
242  }
243 }
244 
245 template <typename V, typename M>
246 void GMRES<V,M>::cycle (V& a_xx, int& a_status, int& a_itcount, RT& a_rnorm0)
247 {
248  BL_PROFILE("GMREA::cycle()");
249 
250  m_res = m_linop->norm2(m_vv[0]);
251  m_grs[0] = m_res;
252 
253  if (m_res == RT(0.0)) {
254  a_status = 0;
255  return;
256  }
257 
258  m_linop->scale(m_vv[0], RT(1.0)/m_res);
259 
260  if (a_itcount == 0) { a_rnorm0 = m_res; }
261 
262  a_status = converged(a_rnorm0,m_res) ? 0 : -1;
263 
264  int it = 0;
265  while (it < m_restrtlen && a_itcount < m_maxiter)
266  {
267  if (m_verbose > 1) {
268  amrex::Print() << "GMRES: iter = " << a_itcount
269  << ", residual = " << m_res << ", " << m_res/a_rnorm0
270  << " (rel.)\n";
271  }
272 
273  if (a_status == 0) { break; }
274 
275  while (m_vv.size() < it+2) {
276  m_vv.emplace_back(m_linop->makeVecRHS());
277  }
278 
279  auto const& vv_it = m_vv[it ];
280  auto & vv_it1 = m_vv[it+1];
281 
282  m_linop->precond(*m_v_tmp_lhs, vv_it);
283  m_linop->apply(vv_it1, *m_v_tmp_lhs);
284 
285  gram_schmidt_orthogonalization(it);
286 
287  auto tt = m_linop->norm2(vv_it1);
288 
289  auto const small = RT((sizeof(RT) == 8) ? 1.e-99 : 1.e-30);
290  bool happyend = (tt < small);
291  if (!happyend) {
292  m_linop->scale(vv_it1, RT(1.0)/tt);
293  }
294 
295  m_hh (it+1,it) = tt;
296  m_hes(it+1,it) = tt;
297 
298  update_hessenberg(it, happyend, m_res);
299 
300  ++it;
301  ++a_itcount;
302  a_status = converged(a_rnorm0, m_res) ? 0 : -1;
303  if (happyend) { break; }
304  }
305 
306  if ((m_verbose > 1) && (a_status != 0 || a_itcount >= m_maxiter)) {
307  amrex::Print() << "GMRES: iter = " << a_itcount
308  << ", residual = " << m_res << ", " << m_res/a_rnorm0
309  << " (rel.)\n";
310  }
311 
312  build_solution(a_xx, it-1);
313 }
314 
315 template <typename V, typename M>
317 {
318  // Two unmodified Gram-Schmidt Orthogonalization
319 
320  BL_PROFILE("GMRES::GramSchmidt");
321 
322  auto& vv_1 = m_vv[it+1];
323 
324  Vector<RT> lhh(it+1);
325 
326  for (int j = 0; j <= it; ++j) {
327  m_hh (j,it) = RT(0.0);
328  m_hes(j,it) = RT(0.0);
329  }
330 
331  for (int ncnt = 0; ncnt < 2 ; ++ncnt)
332  {
333  for (int j = 0; j <= it; ++j) {
334  lhh[j] = m_linop->dotProduct(vv_1, m_vv[j]);
335  }
336 
337  for (int j = 0; j <= it; ++j) {
338  m_linop->increment(vv_1, m_vv[j], -lhh[j]);
339  m_hh (j,it) += lhh[j];
340  m_hes(j,it) -= lhh[j];
341  }
342  }
343 }
344 
345 template <typename V, typename M>
346 void GMRES<V,M>::update_hessenberg (int it, bool happyend, RT& res)
347 {
348  BL_PROFILE("GMRES::update_hessenberg()");
349 
350  for (int j = 1; j <= it; ++j) {
351  auto tt = m_hh(j-1,it);
352  m_hh(j-1,it) = m_cc[j-1] * tt + m_ss[j-1] * m_hh(j,it);
353  m_hh(j ,it) = m_cc[j-1] * m_hh(j,it) - m_ss[j-1] * tt;
354  }
355 
356  if (!happyend)
357  {
358  auto tt = std::sqrt(m_hh(it,it)*m_hh(it,it) + m_hh(it+1,it)*m_hh(it+1,it));
359  m_cc[it] = m_hh(it ,it) / tt;
360  m_ss[it] = m_hh(it+1,it) / tt;
361  m_grs[it+1] = - (m_ss[it] * m_grs[it]);
362  m_grs[it ] = m_cc[it] * m_grs[it];
363  m_hh(it,it) = m_cc[it] * m_hh(it,it) + m_ss[it] * m_hh(it+1,it);
364  res = std::abs(m_grs[it+1]);
365  }
366  else
367  {
368  res = RT(0.0);
369  }
370 }
371 
372 template <typename V, typename M>
373 void GMRES<V,M>::build_solution (V& a_xx, int const it)
374 {
375  BL_PROFILE("GMRES:build_solution()");
376 
377  if (it < 0) { return; }
378 
379  if (m_hh(it,it) != RT(0.0)) {
380  m_grs[it] /= m_hh(it,it);
381  } else {
382  m_grs[it] = RT(0.0);
383  }
384 
385  for (int ii = 1; ii <= it; ++ii) {
386  int k = it - ii;
387  auto tt = m_grs[k];
388  for (int j = k+1; j <= it; ++j) {
389  tt -= m_hh(k,j) * m_grs[j];
390  }
391  m_grs[k] = tt / m_hh(k,k);
392  }
393 
394  m_linop->setToZero(*m_v_tmp_rhs);
395  for (int ii = 0; ii < it+1; ++ii) {
396  m_linop->increment(*m_v_tmp_rhs, m_vv[ii], m_grs[ii]);
397  }
398 
399  m_linop->precond(*m_v_tmp_lhs, *m_v_tmp_rhs);
400  m_linop->increment(a_xx, *m_v_tmp_lhs, RT(1.0));
401 }
402 
403 template <typename V, typename M>
404 void GMRES<V,M>::compute_residual (V& a_rr, V const& a_xx, V const& a_bb)
405 {
406  BL_PROFILE("GMRES::compute_residual()");
407  m_linop->assign(*m_v_tmp_lhs, a_xx);
408  m_linop->apply(*m_v_tmp_rhs, *m_v_tmp_lhs);
409  m_linop->linComb(a_rr, RT(1.0), a_bb, RT(-1.0), *m_v_tmp_rhs);
410 }
411 
412 }
413 #endif
#define BL_PROFILE(a)
Definition: AMReX_BLProfiler.H:551
#define AMREX_ALWAYS_ASSERT(EX)
Definition: AMReX_BLassert.H:50
GMRES.
Definition: AMReX_GMRES.H:61
int getNumIters() const
Gets the number of iterations.
Definition: AMReX_GMRES.H:94
void solve(V &a_sol, V const &a_rhs, RT a_tol_rel, RT a_tol_abs, int a_its=-1)
Solve the linear system.
Definition: AMReX_GMRES.H:193
int m_status
Definition: AMReX_GMRES.H:117
void build_solution(V &a_xx, int it)
Definition: AMReX_GMRES.H:373
void allocate_scratch()
Definition: AMReX_GMRES.H:142
int m_verbose
Definition: AMReX_GMRES.H:114
void define(M &linop)
Definition: AMReX_GMRES.H:168
RT m_rtol
Definition: AMReX_GMRES.H:120
bool converged(RT r0, RT r) const
Definition: AMReX_GMRES.H:187
Vector< RT > m_hh_1d
Definition: AMReX_GMRES.H:122
RT m_atol
Definition: AMReX_GMRES.H:121
int m_maxiter
Definition: AMReX_GMRES.H:115
void cycle(V &a_xx, int &a_status, int &a_itcount, RT &a_rnorm0)
Definition: AMReX_GMRES.H:246
Table2D< RT > m_hes
Definition: AMReX_GMRES.H:125
Table2D< RT > m_hh
Definition: AMReX_GMRES.H:124
void setVerbose(int v)
Sets verbosity.
Definition: AMReX_GMRES.H:85
void compute_residual(V &a_rr, V const &a_xx, V const &a_bb)
Definition: AMReX_GMRES.H:404
Vector< RT > m_cc
Definition: AMReX_GMRES.H:127
M * m_linop
Definition: AMReX_GMRES.H:132
int m_its
Definition: AMReX_GMRES.H:116
int m_restrtlen
Definition: AMReX_GMRES.H:118
void setMaxIters(int niters)
Sets the max number of iterations.
Definition: AMReX_GMRES.H:91
std::unique_ptr< V > m_v_tmp_lhs
Definition: AMReX_GMRES.H:130
Vector< RT > m_hes_1d
Definition: AMReX_GMRES.H:123
void setRestartLength(int rl)
Sets restart length. The default is 30.
Definition: AMReX_GMRES.H:158
void update_hessenberg(int it, bool happyend, RT &res)
Definition: AMReX_GMRES.H:346
void gram_schmidt_orthogonalization(int it)
Definition: AMReX_GMRES.H:316
RT getResidualNorm() const
Gets the 2-norm of the residual.
Definition: AMReX_GMRES.H:100
GMRES()
Definition: AMReX_GMRES.H:136
Vector< V > m_vv
Definition: AMReX_GMRES.H:131
Vector< RT > m_grs
Definition: AMReX_GMRES.H:126
Vector< RT > m_ss
Definition: AMReX_GMRES.H:128
void clear()
Definition: AMReX_GMRES.H:175
int getStatus() const
Gets the solver status.
Definition: AMReX_GMRES.H:97
std::unique_ptr< V > m_v_tmp_rhs
Definition: AMReX_GMRES.H:129
RT m_res
Definition: AMReX_GMRES.H:119
typename M::RT RT
Definition: AMReX_GMRES.H:64
This class provides the user with a few print options.
Definition: AMReX_Print.H:35
@ max
Definition: AMReX_ParallelReduce.H:17
static constexpr int M
Definition: AMReX_OpenBC.H:13
Definition: AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T abs(const GpuComplex< T > &a_z) noexcept
Return the absolute value of a complex number.
Definition: AMReX_GpuComplex.H:356
double second() noexcept
Definition: AMReX_Utility.cpp:922
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
Return the square root of a complex number.
Definition: AMReX_GpuComplex.H:373