Block-Structured AMR Software Framework
AMReX_PCGSolver.H
Go to the documentation of this file.
1 #ifndef AMREX_PCG_SOLVER_H_
2 #define AMREX_PCG_SOLVER_H_
3 #include <AMReX_Config.H>
4 
5 #include <AMReX_Algorithm.H>
6 #include <AMReX_Array.H>
7 #include <cmath>
8 #include <type_traits>
9 
10 namespace amrex {
11 
22 template <int N, typename T, typename M, typename P>
25  M const& mat, P const& precond, int maxiter, T rel_tol)
26 {
27  static_assert(std::is_floating_point_v<T>);
28 
29  T rnorm0 = 0;
30  for (int i = 0; i < N; ++i) {
31  rnorm0 = std::max(rnorm0, std::abs(r[i]));
32  }
33  if (rnorm0 == 0) { return 0; }
34 
35  int iter = 0;
36  T rho_prev = T(1.0); // initialized to quiet gcc warning
37  T p[N];
38  for (iter = 1; iter <= maxiter; ++iter) {
39  T z[N];
40  precond(z, r);
41  T rho = 0;
42  for (int i = 0; i < N; ++i) { rho += r[i]*z[i]; }
43  if (rho == 0) { break; }
44  if (iter == 1) {
45  for (int i = 0; i < N; ++i) { p[i] = z[i]; }
46  } else {
47  auto rr = rho * (T(1.0)/rho_prev);
48  for (int i = 0; i < N; ++i) {
49  p[i] = z[i] + rr * p[i];
50  }
51  }
52  T q[N];
53  mat(q, p);
54  T pq = 0;
55  for (int i = 0; i < N; ++i) { pq += p[i]*q[i]; }
56  if (pq == 0) { break; }
57  T alpha = rho * (T(1.0)/pq);
58  T rnorm = 0;
59  for (int i = 0; i < N; ++i) {
60  x[i] += alpha * p[i];
61  r[i] -= alpha * q[i];
62  rnorm = std::max(rnorm, std::abs(r[i]));
63  }
64  if (rnorm <= rnorm0*rel_tol) { break; }
65  rho_prev = rho;
66  }
67  return iter;
68 }
69 
70 }
71 
72 #endif
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
#define AMREX_RESTRICT
Definition: AMReX_Extension.H:37
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
@ max
Definition: AMReX_ParallelReduce.H:17
static constexpr int M
Definition: AMReX_OpenBC.H:13
static constexpr int P
Definition: AMReX_OpenBC.H:14
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
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int pcg_solve(T *AMREX_RESTRICT x, T *AMREX_RESTRICT r, M const &mat, P const &precond, int maxiter, T rel_tol)
Preconditioned conjugate gradient solver.
Definition: AMReX_PCGSolver.H:24