diff --git a/programs/petsc-gmres/FOM.awk b/programs/petsc-gmres/FOM.awk
new file mode 100644
index 0000000..14ea4cb
--- /dev/null
+++ b/programs/petsc-gmres/FOM.awk
@@ -0,0 +1,31 @@
+#!/usr/bin/awk -f
+# Extract the per-routine detail behind the FOM from a PETSc -log_view file.
+# Column positions are -log_view's fixed layout: $4 time(max), $6 flop(max),
+# $21 total Mflop/s, $22 GPU Mflop/s. $4/$22 read "n/a" unless the run also
+# passed -log_view_gpu_time, so run.sh passes it on every GPU system.
+# Original by A. Suzuki, 25 Aug. 2026.
+BEGIN{
+   timeMatMult = 0.0; flopMatMult = 0.0; totalflopsMatMult = 0.0; gpuflopsMatMult = 0.0;
+   timeKSPSolve = 0.0; flopKSPSolve = 0.0; totalflopsKSPSolve = 0.0; gpuflopsKSPSolve = 0.0;
+   timeSFPack = 0.0;
+   timeSFUnpack = 0.0;
+}
+# space is mandatory to exclude MatMult{Add,Transpose}
+/^MatMult /{
+   timeMatMult = $4; flopMatMult = $6; totalflopsMatMult = $21; gpuflopsMatMult = $22;
+}
+/^KSPSolve/{
+   timeKSPSolve = $4; flopKSPSolve = $6; totalflopsKSPSolve = $21; gpuflopsKSPSolve = $22;
+}
+/^SFPack/{
+   timeSFPack = $4;
+}
+/^SFUnpack/{
+   timeSFUnpack = $4;
+}
+END{
+   printf("#routine     \ttime(sec) \tflop      flop/s(total) flop/s(GPU)\n");
+   printf("MatMult      \t%s \t%s \t%s \t%s\n", timeMatMult, flopMatMult, totalflopsMatMult, gpuflopsMatMult);
+   printf("KSPSolve     \t%s \t%s \t%s \t%s\n", timeKSPSolve, flopKSPSolve, totalflopsKSPSolve, gpuflopsKSPSolve);
+   printf("SFPack/Unpack\t%s\n", timeSFPack + timeSFUnpack);
+}
diff --git a/programs/petsc-gmres/README.md b/programs/petsc-gmres/README.md
new file mode 100644
index 0000000..7e4ed8c
--- /dev/null
+++ b/programs/petsc-gmres/README.md
@@ -0,0 +1,72 @@
+# petsc-gmres
+
+A PETSc KSP (GMRES + GAMG algebraic multigrid preconditioner) benchmark:
+load the `stokes2` sparse matrix (4,260,568 rows, 256,285,536 nonzeros —
+a Stokes flow saddle-point system on an ellipsoid mesh, generated via
+Gmsh/FreeFEM by Atsushi Suzuki), solve `Ax = b` for a known `x`, report
+the relative L2 error and GMRES-iterate wall-time (`FOM: ranks=<N>
+ksp_iter_time_s=<t>`).
+
+The FOM measures the GMRES **iterate only** — `KSPSetUp(ksp)` is called
+explicitly before the timed `KSPSolve`, so GAMG setup (coarsening + PtAP)
+is excluded. The iterate is the MatMult-bound SpMV loop that tracks memory
+bandwidth / Flop/s; setup is communication-bound and is read off
+`-log_view` (`PCSetUp_GAMG`, `MatMult`, ...) instead. (Previously the FOM
+included setup; it was narrowed after `-log_view` showed setup was ~70% of
+the timed `KSPSolve` on a 4-GPU run.)
+
+Source (`src/GMRES-PETSc.cpp`) is vendored directly in this directory
+rather than fetched at build time. PETSc itself is fetched via
+`bk_fetch_source` from its official repo, pinned to `v3.25.2`.
+
+Confirmed working on Rikyu (GB200 NVL4), Fugaku (A64FX), and R-CCS
+Cloud's DGX Spark (`RC_DGXSP`, GB10 Blackwell) — see this app's `build.sh`
+for the per-system recipe.
+
+## Correctness
+
+`relative L2 norm of the error` for `stokes2` lands near **0.23–0.26** at
+any rank count — higher than the ~0.03 typical for SPD matrices like the
+earlier `audikw_1`, because `stokes2` is a Stokes flow saddle-point system
+(velocity + pressure), which is indefinite and harder to solve. GMRES
+hits the default iteration limit rather than converging tightly, but the
+iteration count is stable across rank counts (193–232 on GPU), so the
+timing comparison is fair. A result far outside that band signals a real
+bug, not benchmark noise.
+
+## `-pc_gamg_square_graph 0`
+
+Required on every system, at any rank count above 1 (single-rank runs
+happen not to trigger it, but don't rely on that). Without it, GAMG's
+default aggressive-coarsening graph-squaring step blows up on the matrix's
+connectivity — a `CUSPARSE_STATUS_INSUFFICIENT_RESOURCES` crash on GPU,
+or a genuine multi-GB single-allocation PETSc "Out of memory" abort on
+CPU. Both `build.sh`'s configure line and `run.sh`'s launch command
+already account for everything needed except this flag, which is passed
+explicitly in `run.sh`.
+
+## Staging the data
+
+The matrix is generated by a FreeFEM script (`StokesPETSc0.edp`) that
+solves a Stokes flow on an ellipsoid mesh and writes the resulting matrix
+directly to PETSc binary format via `ObjectView(A, format="binary")` —
+the same format as PETSc's own `MatView`, so no separate conversion step
+(e.g. `mtx2petsc`) is needed. The `MatLoad` call in `GMRES-PETSc.cpp`
+reads it directly.
+
+Pre-staged locations (same convention as this repo's `ffb` and
+`LQCD_dw_solver`, which pre-stage their own — much larger — source
+archives at a fixed per-system path rather than fetching them at
+build/run time):
+
+| system | path |
+|---|---|
+| RIKYU | `/data1/rkp00015/benchkit-data/petsc-gmres/stokes2.dat` |
+| Fugaku | `/vol0500/data/ra250029/benchkit-data/petsc-gmres/stokes2.dat` — `/vol0002` (the volume the first attempt used) turned out to have a 0-byte hard quota limit for every group tried, `/vol0005` didn't; this app's `FJ` queue.csv template already declares `GFSCACHE` for `/vol0002:/vol0003:/vol0004:/vol0005`, so no per-app `-x` handling was needed once the right volume was found. Use `/vol0500`, not `/vol0005`, in an actual path — `/vol0005` is what a filesystem tool reports as the "resolved" location, but it isn't what resolves from a real compute-node job |
+| RC_DGXSP | `/home/users/william.dawson/benchkit-data/petsc-gmres/stokes2.dat` (this system has no separate group-storage tier, so — unlike RIKYU/Fugaku — this lives under `$HOME`) |
+
+To re-stage on a system: copy the `stokes2.dat` file (3,092,468,720 bytes)
+from an existing location above. To regenerate from scratch, run the
+FreeFEM script (`mpirun -np 8 FreeFem++-mpi StokesPETSc0.edp`) which writes
+`A.dat` in PETSc binary format — rename to `stokes2.dat`. Verify: 4,260,568
+× 4,260,568, nnz = 256,285,536.
diff --git a/programs/petsc-gmres/build.sh b/programs/petsc-gmres/build.sh
new file mode 100644
index 0000000..a418889
--- /dev/null
+++ b/programs/petsc-gmres/build.sh
@@ -0,0 +1,131 @@
+#!/bin/bash
+set -euo pipefail
+
+system="$1"
+
+# BenchKit invokes this as `bash programs/petsc-gmres/build.sh <system>`
+# from the repo root, not from inside this directory -- $PWD is the repo
+# root throughout (matching PETSC_DIR/ARTIFACT_DIR below, and bk_fetch_source's
+# own convention). src/GMRES-PETSc.cpp is *this script's own* source, so
+# anchor it to the script's location (APP_DIR) instead of assuming a
+# caller cwd -- found by actually running this through BenchKit's own
+# invocation convention rather than just replicating its commands by hand.
+APP_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+
+PETSC_REPO="https://gitlab.com/petsc/petsc.git"
+PETSC_TAG="v3.25.2"
+PETSC_DIR="${PWD}/petsc"
+ARTIFACT_DIR="${PWD}/artifacts"
+
+source scripts/bk_functions.sh
+
+mkdir -p "${ARTIFACT_DIR}"
+bk_fetch_source "${PETSC_REPO}" "petsc" "${PETSC_TAG}"
+
+case "$system" in
+  RIKYU)
+    module load nvhpc-hpcx/26.3
+    export PETSC_ARCH=arch-rikyu-cuda
+    NVHPC=/shared/software/hpc_sdk/Linux_aarch64/26.3
+    NVPL="${NVHPC}/math_libs/nvpl/lib"
+    (
+      cd "${PETSC_DIR}"
+      # --with-fc=0: nvhpc-hpcx's Fortran wrapper doesn't support an F2018
+      # pointer-initialization feature PETSc's Fortran bindings need; this
+      # app is C++ only anyway. --with-cuda=1: run.sh launches this system
+      # with -mat_type aijcusparse, which a CPU-only PETSc does not
+      # register -- the two must agree. cuda-arch 100 is sm_100 (B200);
+      # hardcoded rather than probed with nvidia-smi (as RC_DGXSP does)
+      # because RIKYU is a cross-mode system whose build runner is a login
+      # node, so a GPU is not guaranteed to be visible at build time.
+      # Serial NVPL BLAS/LAPACK: this app runs flat MPI, no OpenMP.
+      ./configure \
+        --with-cc=mpicc --with-cxx=mpicxx --with-fc=0 \
+        --with-debugging=0 \
+        --with-cuda=1 --with-cuda-dir="${NVHPC}/cuda/13.1" --with-cuda-arch=100 \
+        --with-blaslapack-lib="-L${NVPL} -Wl,-rpath,${NVPL} -lnvpl_lapack_lp64_seq -lnvpl_blas_lp64_seq" \
+        COPTFLAGS='-O3' CXXOPTFLAGS='-O3' CUDAOPTFLAGS='-O3'
+      make PETSC_DIR="${PETSC_DIR}" PETSC_ARCH="${PETSC_ARCH}" -j8 all
+    )
+    mpicxx -O3 -I"${PETSC_DIR}/include" -I"${PETSC_DIR}/${PETSC_ARCH}/include" \
+      -I"${NVHPC}/cuda/13.1/include" \
+      "${APP_DIR}/src/GMRES-PETSc.cpp" -o "${ARTIFACT_DIR}/GMRES-PETSc" \
+      -Xlinker -rpath="${PETSC_DIR}/${PETSC_ARCH}/lib" \
+      -L"${PETSC_DIR}/${PETSC_ARCH}/lib" -lpetsc \
+      -L"${NVHPC}/cuda/13.1/lib64" -Xlinker -rpath="${NVHPC}/cuda/13.1/lib64" -lcudart
+    ;;
+  Fugaku)
+    # The LLVM cross-compiler (mpiclang++) measured faster than the
+    # Fujitsu compiler here, and Fugaku's Spack-provided PETSc hits the
+    # GAMG square-graph crash as a bare SIGKILL rather than a catchable
+    # PETSc error, so build our own.
+    module load lang/tcsds-1.2.43
+    module load LLVM/llvmorg-22.1.0
+    export PETSC_ARCH=arch-fugaku-llvm
+    TCSDS=/opt/FJSVxtclanga/tcsds-latest/lib64
+    ELF=/opt/FJSVxos/devkit/aarch64/rfs/usr/lib64/libelf.so
+    (
+      cd "${PETSC_DIR}"
+      # Login nodes are x86, compute nodes are A64FX/aarch64 -- a genuine
+      # cross-compile, hence --with-batch. BLAS/LAPACK is Fujitsu's own
+      # fjlapacksve (serial), not OpenBLAS -- OpenBLAS pulls in an
+      # unresolvable chain of Fortran-runtime symbols on this toolchain
+      # (RIKEN's own docs: "use of libraries provided with the Fujitsu
+      # compiler with the other compiler environments is not supported",
+      # but they document this exact non-Fujitsu-compiler linkage recipe).
+      ./configure \
+        --with-cc=mpiclang --with-cxx=mpiclang++ --with-fc=0 \
+        --with-debugging=0 --with-batch \
+        --with-blaslapack-lib="${TCSDS}/libfjlapacksve.so ${TCSDS}/libfj90i.so ${TCSDS}/libfj90f.so ${TCSDS}/libfjsrcinfo.so ${TCSDS}/libfjcrt.so ${ELF}" \
+        --PETSC_ARCH="${PETSC_ARCH}"
+      make PETSC_DIR="${PETSC_DIR}" PETSC_ARCH="${PETSC_ARCH}" all
+    )
+    # PETSCLIB is deliberately unquoted below: it's a space-separated list
+    # of separate -L/-Wl,-rpath/-l flags, not one path -- quoting it (an
+    # earlier version of this script did) passes the whole thing as a
+    # single malformed argument to the linker instead of word-splitting
+    # it into individual flags. Found the same way as the APP_DIR fix
+    # above: by actually running this build, not just reading it.
+    PETSCLIB="$(grep '^PETSC_WITH_EXTERNAL_LIB' "${PETSC_DIR}/${PETSC_ARCH}/lib/petsc/conf/petscvariables" | cut -d= -f2-)"
+    # shellcheck disable=SC2086
+    mpiclang++ -O3 -I"${PETSC_DIR}/include" -I"${PETSC_DIR}/${PETSC_ARCH}/include" \
+      "${APP_DIR}/src/GMRES-PETSc.cpp" -o "${ARTIFACT_DIR}/GMRES-PETSc" \
+      $PETSCLIB
+    ;;
+  RC_DGXSP)
+    # GPU build, 1 rank/GPU. -pc_gamg_square_graph 0 (set in run.sh)
+    # avoids a cuSPARSE crash on this matrix's connectivity during GAMG's
+    # aggressive-coarsening graph-squaring step.
+    source /etc/profile.d/modules.sh
+    module load system/ng-dgx nvhpc-hpcx
+    export PETSC_ARCH=arch-dgxsp-cuda
+    CUDA_ARCH=$(nvidia-smi --query-gpu=compute_cap --format=csv,noheader | head -1 | tr -d '.')
+    MATHLIBS_ROOT="$(dirname "$(dirname "$(dirname "$(command -v nvcc)")")")/math_libs"
+    MATHLIBS_CANDIDATES=("${MATHLIBS_ROOT}"/*/lib64)
+    if [ "${#MATHLIBS_CANDIDATES[@]}" -ne 1 ] || [ ! -d "${MATHLIBS_CANDIDATES[0]}" ]; then
+      echo "Expected one NVHPC math_libs/*/lib64 directory under ${MATHLIBS_ROOT}" >&2
+      exit 1
+    fi
+    MATHLIBS="${MATHLIBS_CANDIDATES[0]}"
+    (
+      cd "${PETSC_DIR}"
+      ./configure \
+        --with-cc=mpicc --with-cxx=mpicxx --with-fc=0 \
+        --with-debugging=0 --with-cuda=1 --with-cuda-arch="${CUDA_ARCH}" \
+        LDFLAGS="-L${MATHLIBS} -Wl,-rpath,${MATHLIBS}" \
+        COPTFLAGS='-O3' CXXOPTFLAGS='-O3'
+      make PETSC_DIR="${PETSC_DIR}" PETSC_ARCH="${PETSC_ARCH}" -j8 all
+    )
+    CUDA_INC=$(dirname "$(command -v nvcc)")/../include
+    mpicxx -O3 -I"${PETSC_DIR}/include" -I"${PETSC_DIR}/${PETSC_ARCH}/include" \
+      -I"${CUDA_INC}" \
+      "${APP_DIR}/src/GMRES-PETSc.cpp" -o "${ARTIFACT_DIR}/GMRES-PETSc" \
+      -Xlinker -rpath="${PETSC_DIR}/${PETSC_ARCH}/lib" \
+      -L"${PETSC_DIR}/${PETSC_ARCH}/lib" -lpetsc \
+      -L"${CUDA_INC}/../lib64" -Xlinker -rpath="${CUDA_INC}/../lib64" -lcudart
+    ;;
+  *)
+    echo "Unknown system: $system" >&2
+    exit 1
+    ;;
+esac
diff --git a/programs/petsc-gmres/list.csv b/programs/petsc-gmres/list.csv
new file mode 100644
index 0000000..f3fe3ce
--- /dev/null
+++ b/programs/petsc-gmres/list.csv
@@ -0,0 +1,7 @@
+system,enable,nodes,numproc_node,nthreads,elapse
+RIKYU,yes,1,1,1,0:30:00
+RIKYU,yes,1,2,1,0:30:00
+RIKYU,yes,1,4,1,0:30:00
+RIKYU,yes,2,4,1,0:30:00
+Fugaku,yes,1,48,1,0:15:00
+RC_DGXSP,yes,1,1,1,0:15:00
diff --git a/programs/petsc-gmres/run.sh b/programs/petsc-gmres/run.sh
new file mode 100644
index 0000000..4fa3ee6
--- /dev/null
+++ b/programs/petsc-gmres/run.sh
@@ -0,0 +1,166 @@
+#!/bin/bash
+set -euo pipefail
+
+system="$1"
+nodes="$2"
+numproc_node="$3"
+nthreads="$4"
+n_ranks=$((nodes * numproc_node))
+
+source scripts/bk_functions.sh
+
+# FOM.awk ships next to this script and BenchKit invokes run.sh from the
+# repo root, so anchor it to the script's own location -- same reason
+# build.sh anchors src/GMRES-PETSc.cpp.
+APP_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+
+ARTIFACT="${PWD}/artifacts/GMRES-PETSc"
+RESULTS_DIR="${PWD}/results"
+mkdir -p "${RESULTS_DIR}"
+: > "${RESULTS_DIR}/result"
+
+if [[ ! -x "${ARTIFACT}" ]]; then
+  echo "Required artifact not found or not executable: ${ARTIFACT}" >&2
+  exit 1
+fi
+
+export OMP_NUM_THREADS="${nthreads}"
+
+# The benchmark matrix (stokes2, 4,260,568 x 4,260,568, 256,285,536 nnz --
+# a Stokes flow saddle-point system from an ellipsoid mesh, generated via
+# Gmsh/FreeFEM and written directly to PETSc binary format by FreeFEM's
+# ObjectView, so no mtx2petsc conversion is needed) is pre-staged at a
+# fixed path per system. See README.md for how each copy was produced and
+# how to re-stage it.
+logfile="solve.log"
+# PETSc's own event log, parsed by FOM.awk for the per-routine breakdown.
+logpetsc="petsc.log"
+
+# GMRES restart (Krylov subspace size). PETSc's default of 30 is too small for
+# these Stokes systems: stokes3/stokes4 stagnate at it and never converge.
+# 100 converges stokes2 at every rank count and stokes3 at 2 and 4 GPUs. The
+# operators are singular, so no restart guarantees convergence and larger is
+# not monotonically better -- 80 converged and then broke down on a rerun of
+# the same config.
+KSP_RESTART="${BK_PETSC_GMRES_RESTART:-100}"
+touch .run_marker
+
+case "${system}" in
+  RIKYU)
+    DATA="${BK_PETSC_GMRES_MATRIX:-/data1/rkp00015/benchkit-data/petsc-gmres/stokes2.dat}"
+    module load nvhpc-hpcx/26.3
+    # One MPI rank per GPU: each rank gets a distinct GPU via
+    # CUDA_VISIBLE_DEVICES (set inside mpirun so OMPI_COMM_WORLD_LOCAL_RANK
+    # is available per-rank), and -mat_type aijcusparse puts the matrix
+    # on-device. Without these the solve runs on CPU even though GPUs
+    # are allocated.
+    mpirun -np "${n_ranks}" -N "${numproc_node}" --bind-to core --map-by core \
+      bash -c 'export CUDA_VISIBLE_DEVICES=$OMPI_COMM_WORLD_LOCAL_RANK; exec "$@"' \
+      _ "${ARTIFACT}" -f "${DATA}" -pc_type gamg -pc_gamg_square_graph 0 \
+      -mat_type aijcusparse -matload_block_size 1 \
+      -ksp_gmres_restart "${KSP_RESTART}" \
+      -log_view ":${logpetsc}" -log_view_gpu_time \
+      > "${logfile}" 2>&1 || true
+    ;;
+  Fugaku)
+    # /vol0002 is at quota (0 byte hard limit -- true for every group
+    # tried), but /vol0005 isn't; this repo's FJ queue.csv template
+    # already declares GFSCACHE for both (and /vol0003, /vol0004), so no
+    # extra -x PJM_LLIO_GFSCACHE handling is needed here.
+    # /vol0500 (not /vol0005 -- the "resolved" path fs_mkdir reported when
+    # this was staged) is what actually resolves from a compute-node job;
+    # found by testing the real run.sh in a real job, not by trusting the
+    # canonical-looking path a filesystem tool reported.
+    DATA="${BK_PETSC_GMRES_MATRIX:-/vol0500/data/ra250029/benchkit-data/petsc-gmres/stokes2.dat}"
+    module load lang/tcsds-1.2.43
+    module load LLVM/llvmorg-22.1.0
+    mpiexec -n "${n_ranks}" \
+      "${ARTIFACT}" -f "${DATA}" -pc_type gamg -pc_gamg_square_graph 0 \
+      -matload_block_size 1 \
+      -ksp_gmres_restart "${KSP_RESTART}" \
+      -log_view ":${logpetsc}" \
+      > "${logfile}" 2>&1 || true
+    # Fugaku's PJM mpiexec writes each rank's real stdout/stderr under
+    # ./output.$PJM_JOBID/, ignoring plain shell redirection for the
+    # application's own output -- fall back to searching for it if the
+    # marker wasn't captured above (same pattern as this repo's sbd).
+    if ! grep -q "^FOM: ranks=" "${logfile}" 2>/dev/null; then
+      found=$(find . -maxdepth 5 -type f -newer .run_marker -name 'stdout*' 2>/dev/null | sort | head -n 1)
+      [[ -n "${found}" ]] && logfile="${found}"
+    fi
+    ;;
+  RC_DGXSP)
+    # GPU run (1 rank/GPU) -- see build.sh. This system has no separate
+    # group-storage tier (see README.md), so the data lives under $HOME
+    # like everything else here.
+    DATA="${BK_PETSC_GMRES_MATRIX:-/home/users/william.dawson/benchkit-data/petsc-gmres/stokes2.dat}"
+    source /etc/profile.d/modules.sh
+    module load system/ng-dgx nvhpc-hpcx
+    mpirun -np "${n_ranks}" \
+      bash -c 'export CUDA_VISIBLE_DEVICES=$OMPI_COMM_WORLD_LOCAL_RANK; exec "$@"' \
+      _ "${ARTIFACT}" -f "${DATA}" -pc_type gamg -pc_gamg_square_graph 0 \
+      -mat_type aijcusparse -matload_block_size 1 \
+      -ksp_gmres_restart "${KSP_RESTART}" \
+      -log_view ":${logpetsc}" -log_view_gpu_time \
+      > "${logfile}" 2>&1 || true
+    ;;
+  *)
+    echo "Unknown system: ${system}" >&2
+    exit 1
+    ;;
+esac
+
+if [[ ! -f "${DATA}" ]]; then
+  echo "Benchmark matrix not found at ${DATA} -- see README.md for how to stage it" >&2
+  exit 1
+fi
+
+if ! grep -q "^FOM: ranks=" "${logfile}" 2>/dev/null; then
+  echo "petsc-gmres success marker not found" >&2
+  echo "---- ${logfile} tail ----" >&2
+  tail -n 80 "${logfile}" >&2 || true
+  exit 1
+fi
+
+ksp_iter_time=$(grep "^FOM: ranks=" "${logfile}" | sed -E 's/.*ksp_iter_time_s=([0-9.]+).*/\1/')
+
+# bk_emit_result must come BEFORE the SECTION lines: scripts/result.sh
+# attaches sections to the FOM block currently open and resets them on each
+# FOM line, so sections emitted first are silently dropped from the Result
+# JSON while result.sh still exits 0.
+bk_emit_result \
+  --fom "${ksp_iter_time}" \
+  --fom-unit s \
+  --fom-version ksp_iter_time \
+  --exp "$(basename "${DATA}" .dat)" \
+  --nodes "${nodes}" \
+  --numproc-node "${numproc_node}" \
+  --nthreads "${nthreads}" >> "${RESULTS_DIR}/result"
+
+# Per-routine breakdown behind the FOM: MatMult and KSPSolve time, flop and
+# flop/s, plus the SF pack+unpack time that stands for halo-exchange cost.
+# Kept as a collected artifact and mirrored into SECTION lines, the only
+# structured channel BenchKit has for sub-timings.
+if [[ -f "${logpetsc}" ]]; then
+  awk -f "${APP_DIR}/FOM.awk" "${logpetsc}" > "${RESULTS_DIR}/fom_details.txt"
+  cat "${RESULTS_DIR}/fom_details.txt"
+  # A time reads "n/a" when -log_view could not attribute wall time to the
+  # event. bk_emit_section rejects a non-numeric time, which under set -e
+  # would abort the run and report failure for an otherwise good result --
+  # warn and skip instead.
+  awk -f "${APP_DIR}/FOM.awk" "${logpetsc}" | awk '
+    function num(v) { return v ~ /^[0-9]+([.][0-9]+)?([eE][-+]?[0-9]+)?$/ }
+    $1 == "MatMult"       { print (num($2) ? "MatMult "         $2 : "!MatMult") }
+    $1 == "KSPSolve"      { print (num($2) ? "KSPSolve "        $2 : "!KSPSolve") }
+    $1 == "SFPack/Unpack" { print (num($2) ? "SFPack_SFUnpack " $2 : "!SFPack_SFUnpack") }' |
+  while read -r _sec _time; do
+    if [[ "${_sec}" == !* ]]; then
+      echo "warning: no numeric time for ${_sec#!} in ${logpetsc}, section not emitted" >&2
+      continue
+    fi
+    bk_emit_section "${_sec}" "${_time}" "" "results/fom_details.txt" \
+      >> "${RESULTS_DIR}/result"
+  done
+else
+  echo "warning: ${logpetsc} not produced, skipping FOM detail breakdown" >&2
+fi
diff --git a/programs/petsc-gmres/src/GMRES-PETSc.cpp b/programs/petsc-gmres/src/GMRES-PETSc.cpp
new file mode 100644
index 0000000..531333b
--- /dev/null
+++ b/programs/petsc-gmres/src/GMRES-PETSc.cpp
@@ -0,0 +1,130 @@
+// CSR data transfered to PETSc  Copyright,  Atsushi Suzuki 30 Jul.2025
+//
+// A PETSc KSP (GMRES + GAMG algebraic multigrid preconditioner) benchmark:
+// load a large sparse SPD matrix, solve Ax = b for a known x, report the
+// relative L2 error and solve wall-time.
+//
+// Loads a pre-converted PETSc binary matrix via MatLoad() against
+// PETSC_COMM_WORLD, so rows are properly partitioned across MPI ranks
+// (memory per rank scales as O(1/ranks), as a real distributed benchmark
+// should). The binary is produced once, offline, from the MatrixMarket
+// form of the benchmark matrix (audikw_1, SuiteSparse Matrix Collection,
+// GHS_psdef group) -- see this app's README.md for where each target
+// system's copy lives and how it was produced.
+static char help[] = "load a PETSc binary matrix and call the KSP solver\n";
+
+#include <petscksp.h>
+
+#if defined(PETSC_HAVE_CUDA)
+  #include <cuda_runtime.h>
+#endif
+
+int main(int argc, char **args)
+{
+  Vec         x, b, u;        /* approx solution, RHS, exact solution */
+  Mat         A;              /* linear system matrix */
+  KSP         ksp;            /* linear solver context */
+  PC          pc;             /* PC context */
+  PetscMPIInt size, rank;
+  PetscBool   flg;
+  PetscScalar one = 1.0;
+  PetscReal   e, e0;
+
+  char        fname[1024];
+  PetscViewer viewer;
+
+  PetscFunctionBeginUser;
+#if defined(PETSC_HAVE_CUDA)
+  // Establish this rank's CUDA context BEFORE PetscInitialize calls MPI_Init.
+  // UCX/UCC probe for CUDA support while bootstrapping MPI_COMM_WORLD; with no
+  // context yet, they cache a rendezvous protocol that deadlocks on the first
+  // inter-node GPU-buffer transfer, hanging any 2+ node GPU run in MatMult's
+  // halo exchange with the GPUs idle. run.sh gives each rank one visible GPU,
+  // so device 0 is this rank's GPU.
+  {
+    cudaError_t cerr = cudaSetDevice(0);
+    if (cerr != cudaSuccess) {
+      fprintf(stderr, "warning: cudaSetDevice(0) failed (%s); no CUDA context "
+                      "before MPI_Init -- multi-node GPU runs may hang\n",
+              cudaGetErrorString(cerr));
+    } else if ((cerr = cudaFree(0)) != cudaSuccess) { /* forces context creation */
+      fprintf(stderr, "warning: cudaFree(0) failed (%s); no CUDA context "
+                      "before MPI_Init -- multi-node GPU runs may hang\n",
+              cudaGetErrorString(cerr));
+    }
+  }
+#endif
+  PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
+  PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
+  PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
+  PetscCall(PetscOptionsGetString(NULL, NULL, "-f", fname, sizeof(fname), &flg));
+  if (!flg) {
+    PetscCall(PetscPrintf(PETSC_COMM_WORLD,
+			  "usage: %s -f <matrix.petscbin> [-ksp_monitor] [-pc_type gamg] [-log_view]\n",
+			  args[0]));
+    PetscCall(PetscFinalize());
+    exit(-1);
+  }
+
+  PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
+  PetscCall(MatSetType(A, MATAIJ));
+  PetscCall(MatSetFromOptions(A));
+  PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, fname, FILE_MODE_READ, &viewer));
+  PetscCall(MatLoad(A, viewer));
+  PetscCall(PetscViewerDestroy(&viewer));
+
+  //     Create and set vectors
+  PetscCall(MatCreateVecs(A, &u, &b));  // create vectors right u and left b
+  PetscCall(VecDuplicate(u, &x));       // copy vector x from u
+  PetscCall(VecSet(u, one));            // all entries of u are one
+  PetscCall(MatMult(A, u, b));          // RHS from b = A * u
+
+  //     Create linear solver context
+  PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
+  PetscCall(KSPSetOperators(ksp, A, A));
+
+  //     Set defulat preconditioner as diagonal preconditioning
+  PetscCall(KSPGetPC(ksp, &pc));
+  PetscCall(PCSetType(pc, PCJACOBI));
+  //     receive options for KSP solver from command line
+  PetscCall(KSPSetFromOptions(ksp));
+
+  // Build the GAMG hierarchy (coarsening + Galerkin/PtAP) explicitly,
+  // OUTSIDE the timed region. The figure of merit is the GMRES iterate --
+  // the MatMult-bound SpMV loop whose cost tracks memory bandwidth / Flop/s
+  // -- not one-shot preconditioner setup, which is communication-bound and
+  // does not characterize the iterative kernel. KSPSolve would otherwise
+  // call KSPSetUp lazily on first use, folding setup into the timed region.
+  // Use -log_view for the per-event breakdown (PCSetUp_GAMG, MatMult, ...).
+  PetscCall(KSPSetUp(ksp));
+
+  PetscLogDouble iter_t0, iter_t1;
+  PetscCall(PetscBarrier((PetscObject)ksp));
+  PetscCall(PetscTime(&iter_t0));
+  PetscCall(KSPSolve(ksp, b, x));
+  PetscCall(PetscBarrier((PetscObject)ksp));
+  PetscCall(PetscTime(&iter_t1));
+  {
+    PetscReal local_dt = (PetscReal)(iter_t1 - iter_t0), max_dt;
+    PetscCallMPI(MPI_Allreduce(&local_dt, &max_dt, 1, MPIU_REAL, MPI_MAX, PETSC_COMM_WORLD));
+    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "FOM: ranks=%d ksp_iter_time_s=%.6f\n", (int)size, (double)max_dt));
+  }
+
+  PetscCall(VecNorm(u, NORM_2, &e0));
+  PetscCall(VecAXPY(x, -1.0, u));
+  PetscCall(VecNorm(x, NORM_2, &e));
+  PetscCall(PetscPrintf(PETSC_COMM_WORLD,
+			"relative L2 norm of the error: %g\n", (double)(e/e0)));
+
+  flg = PETSC_FALSE;
+  PetscCall(PetscOptionsGetBool(NULL, NULL, "-print_error", &flg, NULL));
+  if (flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Infinity norm of the error: %g\n", (double)e));
+
+  PetscCall(KSPDestroy(&ksp));
+  PetscCall(VecDestroy(&u));
+  PetscCall(VecDestroy(&x));
+  PetscCall(VecDestroy(&b));
+  PetscCall(MatDestroy(&A));
+  PetscCall(PetscFinalize());
+  return 0;
+}