PoPS/doxygen/generic__krylov_8hpp_source.html

#pragma once


#include <pops/core/foundation/types.hpp>

#include <pops/mesh/storage/mf_arith.hpp>

#include <pops/mesh/storage/multifab.hpp>

#include <pops/numerics/elliptic/linear/krylov_result.hpp>  // pops::KrylovResult (shared)


#include <cmath>

#include <functional>

#include <stdexcept>

#include <vector>


namespace pops {


// KrylovResult is shared via krylov_result.hpp (included above), so this header and the

// GeometricMG-coupled krylov_solver.hpp never carry hand-synchronised copies of the struct.


using ApplyFn = std::function<void(MultiFab& out, const MultiFab& in)>;


namespace detail {


inline Real krylov_dot(const MultiFab& x, const MultiFab& y) {

  return x.ncomp() > 1 ? dot_all(x, y) : dot(x, y);

}


inline Real krylov_l2_norm(const MultiFab& x) {

  return std::sqrt(krylov_dot(x, x));

}


inline void require_max_iter(int max_iters) {

  if (max_iters <= 0)

    throw std::invalid_argument("dynamic solver loops require max_iter");

}


inline constexpr Real kKrylovTiny = Real(1e-300);


}  // namespace detail


inline KrylovResult richardson_solve(const ApplyFn& A, MultiFab& phi, const MultiFab& rhs,

                                     Real omega, Real rel_tol, int max_iters) {

  detail::require_max_iter(max_iters);

  // Scratch allocated ONCE, co-distributed with phi; reused every iteration (no in-loop alloc).

  MultiFab r(phi.box_array(), phi.dmap(), phi.ncomp(), phi.n_grow());


  const Real bnorm = detail::krylov_l2_norm(rhs);        // COLLECTIVE

  const Real norm0 = bnorm > Real(0) ? bnorm : Real(1);  // zero rhs -> absolute base


  A(r, phi);                               // r = A(phi)

  lincomb(r, Real(1), rhs, Real(-1), r);   // r = b - A(phi)

  Real rnorm = detail::krylov_l2_norm(r);  // COLLECTIVE

  KrylovResult res;

  res.rel_residual = rnorm / norm0;

  if (rnorm <= rel_tol * norm0) {  // already converged (warm start)

    res.converged = true;

    return res;

  }


  for (int k = 1; k <= max_iters; ++k) {

    saxpy(phi, omega, r);                   // phi <- phi + omega r

    A(r, phi);                              // r = A(phi)

    lincomb(r, Real(1), rhs, Real(-1), r);  // r = b - A(phi)

    rnorm = detail::krylov_l2_norm(r);      // COLLECTIVE (all ranks, same value)

    res.iters = k;

    res.rel_residual = rnorm / norm0;

    if (rnorm <= rel_tol * norm0) {

      res.converged = true;

      return res;

    }

  }

  return res;  // max_iters reached: best effort (converged=false)

}


inline KrylovResult cg_solve(const ApplyFn& A, MultiFab& phi, const MultiFab& rhs, Real rel_tol,

                             int max_iters) {

  detail::require_max_iter(max_iters);

  // Scratch allocated ONCE, co-distributed with phi; reused every iteration. p carries phi's ghost

  // width because A(Ap, p) may read p's ghosts; r and Ap need none but share the layout for clarity.

  MultiFab r(phi.box_array(), phi.dmap(), phi.ncomp(), phi.n_grow());

  MultiFab p(phi.box_array(), phi.dmap(), phi.ncomp(), phi.n_grow());

  MultiFab Ap(phi.box_array(), phi.dmap(), phi.ncomp(), phi.n_grow());


  const Real bnorm = detail::krylov_l2_norm(rhs);  // COLLECTIVE

  const Real norm0 = bnorm > Real(0) ? bnorm : Real(1);


  A(r, phi);                               // r = A(phi)

  lincomb(r, Real(1), rhs, Real(-1), r);   // r = b - A(phi)

  Real rnorm = detail::krylov_l2_norm(r);  // COLLECTIVE

  KrylovResult res;

  res.rel_residual = rnorm / norm0;

  if (rnorm <= rel_tol * norm0) {  // already converged (warm start)

    res.converged = true;

    return res;

  }


  lincomb(p, Real(1), r, Real(0), r);      // p_0 = r_0 (copy via lincomb)

  Real rs_old = detail::krylov_dot(r, r);  // COLLECTIVE: ||r||^2 (all components if ncomp>1)


  for (int k = 1; k <= max_iters; ++k) {

    A(Ap, p);                                    // Ap = A(p)

    const Real pAp = detail::krylov_dot(p, Ap);  // COLLECTIVE (all components if ncomp>1)

    if (std::fabs(pAp) < detail::kKrylovTiny) {  // breakdown (A not SPD / null direction)

      res.iters = k - 1;

      res.rel_residual = rnorm / norm0;

      return res;

    }

    const Real alpha = rs_old / pAp;

    saxpy(phi, alpha, p);                          // phi <- phi + alpha p

    saxpy(r, -alpha, Ap);                          // r   <- r - alpha Ap

    const Real rs_new = detail::krylov_dot(r, r);  // COLLECTIVE: new ||r||^2 (all components)

    rnorm = std::sqrt(rs_new);

    res.iters = k;

    res.rel_residual = rnorm / norm0;

    if (rnorm <= rel_tol * norm0) {

      res.converged = true;

      return res;

    }

    const Real beta = rs_new / rs_old;

    lincomb(p, Real(1), r, beta, p);  // p <- r + beta p

    rs_old = rs_new;

  }

  return res;  // max_iters reached: best effort (converged=false)

}


inline KrylovResult bicgstab_solve(const ApplyFn& A, const ApplyFn& precond, MultiFab& phi,

                                   const MultiFab& rhs, Real rel_tol, int max_iters) {

  detail::require_max_iter(max_iters);

  const bool has_precond = static_cast<bool>(precond);


  // Scratch allocated ONCE, co-distributed with phi; reused every iteration (mirror of

  // TensorKrylovSolver's fixed fields). The matvec inputs (p / s, or phat / shat when preconditioned)

  // carry phi's ghost width since A may read their ghosts; the pointwise scratch (r, rhat, v, t)

  // shares the layout. phat/shat are only allocated when a preconditioner is supplied.

  const BoxArray& ba = phi.box_array();

  const DistributionMapping& dm = phi.dmap();

  const int nc = phi.ncomp();

  const int ng = phi.n_grow();

  MultiFab r(ba, dm, nc, ng), rhat(ba, dm, nc, ng), p(ba, dm, nc, ng);

  MultiFab v(ba, dm, nc, ng), s(ba, dm, nc, ng), t(ba, dm, nc, ng);

  MultiFab phat, shat;

  if (has_precond) {

    phat = MultiFab(ba, dm, nc, ng);

    shat = MultiFab(ba, dm, nc, ng);

  }


  const Real bnorm = detail::krylov_l2_norm(rhs);  // COLLECTIVE

  const Real norm0 = bnorm > Real(0) ? bnorm : Real(1);


  A(v, phi);                               // v = A(phi) (scratch reuse for r0)

  lincomb(r, Real(1), rhs, Real(-1), v);   // r = b - A(phi)

  Real rnorm = detail::krylov_l2_norm(r);  // COLLECTIVE

  KrylovResult res;

  res.rel_residual = rnorm / norm0;

  if (rnorm <= rel_tol * norm0) {  // already converged (warm start)

    res.converged = true;

    return res;

  }


  lincomb(rhat, Real(1), r, Real(0), r);  // rhat = frozen r0 (shadow vector)

  p.set_val(Real(0));

  v.set_val(Real(0));

  Real rho_prev = Real(1), alpha = Real(1), omega = Real(1);


  for (int k = 1; k <= max_iters; ++k) {

    const Real rho = detail::krylov_dot(rhat, r);  // COLLECTIVE (all components if ncomp>1)

    if (std::fabs(rho) < detail::kKrylovTiny || std::fabs(omega) < detail::kKrylovTiny) {

      res.iters = k - 1;  // breakdown: best effort

      res.rel_residual = rnorm / norm0;

      return res;

    }

    const Real beta = (rho / rho_prev) * (alpha / omega);

    lincomb(p, Real(1), p, -omega, v);            // p <- p - omega v

    lincomb(p, beta, p, Real(1), r);              // p <- r + beta p

    MultiFab& phat_ref = has_precond ? phat : p;  // identity precond: phat = p

    if (has_precond)

      precond(phat, p);                                   // phat = M^{-1} p

    A(v, phat_ref);                                       // v = A(phat)

    const Real rhat_dot_v = detail::krylov_dot(rhat, v);  // COLLECTIVE (all components if ncomp>1)

    if (std::fabs(rhat_dot_v) < detail::kKrylovTiny) {

      res.iters = k - 1;

      res.rel_residual = rnorm / norm0;

      return res;

    }

    alpha = rho / rhat_dot_v;

    lincomb(s, Real(1), r, -alpha, v);             // s <- r - alpha v

    saxpy(phi, alpha, phat_ref);                   // phi <- phi + alpha phat (partial)

    const Real snorm = detail::krylov_l2_norm(s);  // COLLECTIVE

    if (snorm <= rel_tol * norm0) {                // mid-iteration convergence

      rnorm = snorm;

      res.iters = k;

      res.rel_residual = rnorm / norm0;

      res.converged = true;

      return res;

    }

    MultiFab& shat_ref = has_precond ? shat : s;  // identity precond: shat = s

    if (has_precond)

      precond(shat, s);                        // shat = M^{-1} s

    A(t, shat_ref);                            // t = A(shat)

    const Real tt = detail::krylov_dot(t, t);  // COLLECTIVE (all components if ncomp>1)

    omega = tt > detail::kKrylovTiny ? detail::krylov_dot(t, s) / tt : Real(0);  // COLLECTIVE

    saxpy(phi, omega, shat_ref);        // phi <- phi + omega shat

    lincomb(r, Real(1), s, -omega, t);  // r <- s - omega t

    rnorm = detail::krylov_l2_norm(r);  // COLLECTIVE

    res.iters = k;

    res.rel_residual = rnorm / norm0;

    if (rnorm <= rel_tol * norm0) {

      res.converged = true;

      return res;

    }

    rho_prev = rho;

  }

  return res;  // max_iters reached: best effort (converged=false)

}


inline KrylovResult gmres_solve(const ApplyFn& A, const ApplyFn& precond, MultiFab& phi,

                                const MultiFab& rhs, Real rel_tol, int max_iters,

                                int restart = 30) {

  detail::require_max_iter(max_iters);

  constexpr int kGmresRestartMax = 50;  // caps the stack Hessenberg system (restart+1 entries)

  const int m = restart < 1 ? 1 : (restart > kGmresRestartMax ? kGmresRestartMax : restart);

  const bool has_precond = static_cast<bool>(precond);


  const BoxArray& ba = phi.box_array();

  const DistributionMapping& dm = phi.dmap();

  const int nc = phi.ncomp();

  const int ng = phi.n_grow();


  // Scratch allocated ONCE, co-distributed with phi; reused across every restart cycle. The Arnoldi

  // basis V holds m+1 vectors; w is the matvec / MGS workspace (carries phi's ghost width since A reads

  // its ghosts); r is the (preconditioned) residual; Mb is M^{-1} b for the relative-residual base.

  std::vector<MultiFab> V;

  V.reserve(static_cast<std::size_t>(m) + 1);

  for (int i = 0; i <= m; ++i)

    V.emplace_back(ba, dm, nc, ng);

  MultiFab w(ba, dm, nc, ng);

  MultiFab r(ba, dm, nc, ng);

  MultiFab Mb;

  if (has_precond)

    Mb = MultiFab(ba, dm, nc, ng);


  // Preconditioned right-hand-side norm: the left-preconditioned residual is M^{-1}(b - A x), so its

  // base is ||M^{-1} b||. With the identity Mb aliases rhs and this is the true ||b||, matching CG.

  if (has_precond)

    precond(Mb, rhs);

  const MultiFab& Mb_ref = has_precond ? Mb : rhs;

  const Real bnorm = detail::krylov_l2_norm(Mb_ref);  // COLLECTIVE

  const Real norm0 = bnorm > Real(0) ? bnorm : Real(1);


  // Fixed-size stack Hessenberg least-squares: H column-major (m columns of m+1 rows), the Givens

  // rotation cosines/sines, the rotated rhs g (beta e1 after rotations), and the back-substituted y.

  // restart <= kGmresRestartMax keeps every array bounded; no heap, no Eigen / LAPACK.

  Real H[kGmresRestartMax + 1][kGmresRestartMax];

  Real cs[kGmresRestartMax];

  Real sn[kGmresRestartMax];

  Real g[kGmresRestartMax + 1];

  Real y[kGmresRestartMax];


  KrylovResult res;

  int total_iters = 0;


  while (total_iters < max_iters) {

    // r = M^{-1}(b - A phi): residual of the current iterate, then left-preconditioned.

    A(w, phi);                              // w = A(phi)

    lincomb(w, Real(1), rhs, Real(-1), w);  // w = b - A(phi)

    if (has_precond)

      precond(r, w);  // r = M^{-1}(b - A phi)

    else

      lincomb(r, Real(1), w, Real(0), w);   // r = b - A phi (copy via lincomb)

    Real beta = detail::krylov_l2_norm(r);  // COLLECTIVE

    res.rel_residual = beta / norm0;

    if (beta <= rel_tol * norm0) {  // converged (also the warm-start early-out)

      res.converged = true;

      return res;

    }


    lincomb(V[0], Real(1) / beta, r, Real(0), r);  // v_0 = r / ||r||

    g[0] = beta;

    for (int i = 1; i <= m; ++i)

      g[i] = Real(0);


    int k = 0;  // Krylov steps taken this cycle

    for (; k < m && total_iters < max_iters; ++k) {

      ++total_iters;

      // Arnoldi: w = M^{-1} A v_k, orthogonalised against v_0..v_k by MODIFIED Gram-Schmidt.

      if (has_precond) {

        A(r, V[k]);     // r = A v_k (reuse r as the unpreconditioned matvec scratch)

        precond(w, r);  // w = M^{-1} A v_k

      } else {

        A(w, V[k]);  // w = A v_k

      }

      for (int i = 0; i <= k; ++i) {

        H[i][k] = detail::krylov_dot(w, V[i]);  // COLLECTIVE (all components if ncomp>1)

        saxpy(w, -H[i][k], V[i]);               // w <- w - H[i][k] v_i

      }

      H[k + 1][k] = detail::krylov_l2_norm(w);  // COLLECTIVE

      if (H[k + 1][k] > detail::kKrylovTiny)

        lincomb(V[k + 1], Real(1) / H[k + 1][k], w, Real(0), w);  // v_{k+1} = w / h

      // else: lucky breakdown (the exact solution lies in the current subspace); v_{k+1} stays unused.


      // Apply the previous Givens rotations to the new Hessenberg column, then a fresh rotation to

      // annihilate H[k+1][k]; carry the rotation through g (the rotated beta e1).

      for (int i = 0; i < k; ++i) {

        const Real t = cs[i] * H[i][k] + sn[i] * H[i + 1][k];

        H[i + 1][k] = -sn[i] * H[i][k] + cs[i] * H[i + 1][k];

        H[i][k] = t;

      }

      const Real denom = std::sqrt(H[k][k] * H[k][k] + H[k + 1][k] * H[k + 1][k]);

      if (denom > detail::kKrylovTiny) {

        cs[k] = H[k][k] / denom;

        sn[k] = H[k + 1][k] / denom;

      } else {  // degenerate column: identity rotation (no division by ~0)

        cs[k] = Real(1);

        sn[k] = Real(0);

      }

      H[k][k] = cs[k] * H[k][k] + sn[k] * H[k + 1][k];

      H[k + 1][k] = Real(0);

      const Real g_next = -sn[k] * g[k];

      g[k] = cs[k] * g[k];

      g[k + 1] = g_next;


      const Real resid = std::fabs(g[k + 1]);  // ||M^{-1}(b - A x)|| WITHOUT a matvec

      res.iters = total_iters;

      res.rel_residual = resid / norm0;

      if (resid <= rel_tol * norm0) {

        ++k;  // include this step in the back-substitution

        break;

      }

    }


    // Back-substitution: solve the k x k upper-triangular R y = g, then phi <- phi + sum y_i v_i.

    for (int i = k - 1; i >= 0; --i) {

      Real sum = g[i];

      for (int j = i + 1; j < k; ++j)

        sum -= H[i][j] * y[j];

      y[i] = std::fabs(H[i][i]) > detail::kKrylovTiny ? sum / H[i][i] : Real(0);

    }

    for (int i = 0; i < k; ++i)

      saxpy(phi, y[i], V[i]);  // phi <- phi + y_i v_i


    // True residual of the updated iterate decides convergence / continuation (the rotated estimate

    // g[k] guides the inner cycle; this recompute keeps the stopping test honest across restarts).

    A(w, phi);                              // w = A(phi)

    lincomb(w, Real(1), rhs, Real(-1), w);  // w = b - A(phi)

    if (has_precond)

      precond(r, w);

    else

      lincomb(r, Real(1), w, Real(0), w);

    const Real rnorm = detail::krylov_l2_norm(r);  // COLLECTIVE

    res.rel_residual = rnorm / norm0;

    if (rnorm <= rel_tol * norm0) {

      res.converged = true;

      return res;

    }

  }

  return res;  // max_iters reached: best effort (converged=false)

}


}  // namespace pops

pops::BoxArray
Ordered list of boxes tiling a level.
Definition box_array.hpp:22

pops::DistributionMapping
Owning MPI rank of each box, indexed by GLOBAL box index (parallel to a BoxArray).
Definition distribution_mapping.hpp:19

pops::MultiFab
Field distributed over a level: decomposition (BoxArray) + distribution (DistributionMapping) + ncomp...
Definition multifab.hpp:33

pops::MultiFab::ncomp
int ncomp() const
Number of components.
Definition multifab.hpp:60

pops::MultiFab::dmap
const DistributionMapping & dmap() const
GLOBAL distribution (owner rank per box).
Definition multifab.hpp:58

pops::MultiFab::n_grow
int n_grow() const
Number of ghost layers.
Definition multifab.hpp:62

pops::MultiFab::box_array
const BoxArray & box_array() const
GLOBAL decomposition of the level (all boxes, all ranks).
Definition multifab.hpp:56

pops::MultiFab::set_val
void set_val(Real v)
Fills all cells (valid + ghosts) of every local fab with v.
Definition multifab.hpp:85

krylov_result.hpp
KrylovResult – the shared result type of every Krylov solve in include/pops/numerics/elliptic/linear.

mf_arith.hpp
MultiFab arithmetic (saxpy, lincomb, norm_inf, dot) over VALID cells.

multifab.hpp
MultiFab: a field DISTRIBUTED over a level (equivalent of AMReX's MultiFab).

pops::detail::kKrylovTiny
constexpr Real kKrylovTiny
Tiny breakdown guard for the BiCGStab scalar recurrences (division by ~0).
Definition generic_krylov.hpp:93

pops::detail::krylov_l2_norm
Real krylov_l2_norm(const MultiFab &x)
GLOBAL L2 norm sqrt(sum x.x), collective (all_reduce_sum).
Definition generic_krylov.hpp:82

pops::detail::krylov_dot
Real krylov_dot(const MultiFab &x, const MultiFab &y)
Krylov inner product x.y, COLLECTIVE (all_reduce_sum).
Definition generic_krylov.hpp:76

pops::detail::require_max_iter
void require_max_iter(int max_iters)
Guards a dynamic solver loop: a non-positive iteration budget is a configuration error.
Definition generic_krylov.hpp:87

pops
Definition amr_hierarchy.hpp:29

pops::gmres_solve
KrylovResult gmres_solve(const ApplyFn &A, const ApplyFn &precond, MultiFab &phi, const MultiFab &rhs, Real rel_tol, int max_iters, int restart=30)
Left-preconditioned restarted GMRES(m), solving A x = b for a GENERAL (possibly NON-symmetric) operat...
Definition generic_krylov.hpp:344

pops::dot_all
Real dot_all(const MultiFab &x, const MultiFab &y)
FULL-component dot Sum_{cells, c} x(.,.,c) * y(.,.,c) over ALL components, reduced over ALL ranks (al...
Definition mf_arith.hpp:182

pops::Real
double Real
Definition types.hpp:30

pops::cg_solve
KrylovResult cg_solve(const ApplyFn &A, MultiFab &phi, const MultiFab &rhs, Real rel_tol, int max_iters)
Conjugate Gradient, solving A x = b for an SPD operator A.
Definition generic_krylov.hpp:155

pops::dot
Real dot(const MultiFab &x, const MultiFab &y, int comp=0)
Dot product Sum_cells x.y over component comp, reduced over ALL ranks (all_reduce).
Definition mf_arith.hpp:163

pops::richardson_solve
KrylovResult richardson_solve(const ApplyFn &A, MultiFab &phi, const MultiFab &rhs, Real omega, Real rel_tol, int max_iters)
Richardson iteration x <- x + omega (b - A x), solving A x = b.
Definition generic_krylov.hpp:108

pops::saxpy
void saxpy(MultiFab &y, Real a, const MultiFab &x)
y <- y + a x over ALL components of the valid cells. Identical layouts required.
Definition mf_arith.hpp:102

pops::sum
Real sum(const MultiFab &mf, int comp=0)
Sum of the VALID cells of component comp, reduced over ALL ranks (all_reduce).
Definition multifab.hpp:121

pops::lincomb
void lincomb(MultiFab &z, Real a, const MultiFab &x, Real b, const MultiFab &y)
z <- a x + b y over ALL components of the valid cells. Identical layouts; aliasing safe.
Definition mf_arith.hpp:133

pops::ApplyFn
std::function< void(MultiFab &out, const MultiFab &in)> ApplyFn
Matrix-free operator callback: out <- A(in).
Definition generic_krylov.hpp:67

pops::bicgstab_solve
KrylovResult bicgstab_solve(const ApplyFn &A, const ApplyFn &precond, MultiFab &phi, const MultiFab &rhs, Real rel_tol, int max_iters)
Preconditioned BiCGStab, solving A x = b for a general (possibly non-symmetric) operator A.
Definition generic_krylov.hpp:223

pops::KrylovResult
Outcome of a Krylov solve: iterations performed, final relative residual, convergence flag.
Definition krylov_result.hpp:16

pops::KrylovResult::iters
int iters
number of iterations performed
Definition krylov_result.hpp:17

pops::KrylovResult::rel_residual
Real rel_residual
||r_final|| / ||b|| (global L2 norm; base 1 when ||b|| == 0)
Definition krylov_result.hpp:18

pops::KrylovResult::converged
bool converged
true if ||r|| <= rel_tol * ||b|| was reached
Definition krylov_result.hpp:19

types.hpp
Base scalar types and the POPS_HD macro (host+device portability).