diff --git a/programs/petsc-gmres/README.md b/programs/petsc-gmres/README.md
new file mode 100644
index 0000000..71793e3
--- /dev/null
+++ b/programs/petsc-gmres/README.md
@@ -0,0 +1,74 @@
+# petsc-gmres
+
+A PETSc KSP (GMRES + GAMG algebraic multigrid preconditioner) benchmark:
+load the `audikw_1` sparse SPD matrix ([SuiteSparse Matrix
+Collection](https://sparse.tamu.edu/GHS_psdef/audikw_1), `GHS_psdef`
+group — a real structural-engineering FEM problem, 943,695 rows,
+77,651,847 nonzeros), solve `Ax = b` for a known `x`, report the relative
+L2 error and solve wall-time (`FOM: ranks=<N> solve_time_s=<t>`).
+
+Source (`src/GMRES-PETSc.cpp`) is vendored directly in this directory
+rather than fetched from a separate repo at build time, since its
+upstream development happens in a private RIKEN-RCCS repository. PETSc
+itself is fetched normally 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. Full strong-scaling results, root-cause
+diagnoses for a couple of real bugs found getting each port working
+(a matrix-distribution bug causing OOM, a GAMG coarsening crash on this
+matrix's connectivity, an MPI-transport gotcha on one system), and the
+underlying build recipes are documented in more depth in that same
+internal repository — not linked here since it isn't publicly readable,
+but available to RIKEN-RCCS members on request.
+
+## Correctness
+
+`relative L2 norm of the error` should land near 0.03–0.04 at any rank
+count (GMRES's default relative-residual tolerance, not a tight solve —
+this is a benchmark, not a production accuracy target). A result 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
+`audikw_1`'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 pre-converted once (offline, not part of `build.sh`/`run.sh`)
+from MatrixMarket format to PETSc's binary format via a small one-time
+conversion tool (also part of the private upstream repo, not shipped
+here — it's not needed at benchmark build/run time, only to produce the
+staged file below once). This avoids every rank re-parsing a
+multi-hundred-MB text file at load time, which was the actual root cause
+of the original reported "runs out of memory for no reason" bug this
+benchmark exists to catch a regression of.
+
+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/audikw_1.petscbin` |
+| Fugaku | `/vol0500/data/ra250029/benchkit-data/petsc-gmres/audikw_1.petscbin` — `/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/audikw_1.petscbin` (this system has no separate group-storage tier, so — unlike RIKYU/Fugaku — this lives under `$HOME`) |
+
+To re-stage on a system with an existing PETSc install: download
+`audikw_1.mtx` from the SuiteSparse Matrix Collection link above, then use
+PETSc's own `MatLoad`/`MatView` binary-viewer round trip (or the private
+repo's `mtx2petsc` tool, if you have access) to write it to PETSc binary
+format. Verify the result: 943,695 × 943,695, nnz = 77,651,847 — the
+standard SuiteSparse download is MatrixMarket `symmetric` format (one
+triangle + diagonal only), so a correct converter must mirror off-diagonal
+entries; a naive read of the file as-is will silently produce a wrong,
+singular matrix with no error.
diff --git a/programs/petsc-gmres/build.sh b/programs/petsc-gmres/build.sh
new file mode 100644
index 0000000..6573a35
--- /dev/null
+++ b/programs/petsc-gmres/build.sh
@@ -0,0 +1,119 @@
+#!/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
+    NVPL=/shared/software/hpc_sdk/Linux_aarch64/26.3/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=0: this is the CPU-only build
+      # (see petsc-benchmarking in RIKEN-RCCS/block1-eeas for why a
+      # +cuda-enabled build's PetscInitialize scales ~linearly with rank
+      # count from CUDA-library dynamic-linking contention, unrelated to
+      # solve performance). Serial NVPL BLAS/LAPACK: this app runs flat
+      # MPI, no OpenMP threading.
+      ./configure \
+        --with-cc=mpicc --with-cxx=mpicxx --with-fc=0 \
+        --with-debugging=0 --with-cuda=0 \
+        --with-blaslapack-lib="-L${NVPL} -Wl,-rpath,${NVPL} -lnvpl_lapack_lp64_seq -lnvpl_blas_lp64_seq" \
+        COPTFLAGS='-O3' CXXOPTFLAGS='-O3'
+      make PETSC_DIR="${PETSC_DIR}" PETSC_ARCH="${PETSC_ARCH}" -j8 all
+    )
+    mpicxx -O3 -I"${PETSC_DIR}/include" -I"${PETSC_DIR}/${PETSC_ARCH}/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
+    ;;
+  Fugaku)
+    # See petsc-build in RIKEN-RCCS/block1-eeas: LLVM cross-compiler
+    # (mpiclang++) beats the Fujitsu compiler here, and Fugaku's own
+    # Spack-provided PETSc hits the same GAMG square-graph crash as a bare
+    # SIGKILL instead of a catchable PETSc error, so we 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 -- audikw_1 is solved on-GPU here (1 rank/GPU), matching
+    # the cross-machine reproduction in RIKEN-RCCS/block1-eeas's
+    # petsc-build. -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=$(dirname "$(command -v nvcc)")/../../math_libs/*/lib64
+    (
+      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
+    )
+    mpicxx -O3 -I"${PETSC_DIR}/include" -I"${PETSC_DIR}/${PETSC_ARCH}/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
+    ;;
+  *)
+    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..5b9095c
--- /dev/null
+++ b/programs/petsc-gmres/list.csv
@@ -0,0 +1,4 @@
+system,enable,nodes,numproc_node,nthreads,elapse
+RIKYU,yes,1,4,1,0:15: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..31fab96
--- /dev/null
+++ b/programs/petsc-gmres/run.sh
@@ -0,0 +1,104 @@
+#!/bin/bash
+set -euo pipefail
+
+system="$1"
+nodes="$2"
+numproc_node="$3"
+nthreads="$4"
+n_ranks=$((nodes * numproc_node))
+
+source scripts/bk_functions.sh
+
+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 (audikw_1, SuiteSparse GHS_psdef group, 943,695 x
+# 943,695, 77,651,847 nnz, converted once to PETSc binary format) is
+# pre-staged at a fixed path per system rather than fetched at build/run
+# time -- same convention as this repo's ffb and LQCD_dw_solver, which
+# pre-stage their (much larger) source archives the same way. See
+# README.md for exactly how each copy was produced and how to re-stage it.
+logfile="solve.log"
+touch .run_marker
+
+case "${system}" in
+  RIKYU)
+    DATA=/data1/rkp00015/benchkit-data/petsc-gmres/audikw_1.petscbin
+    module load nvhpc-hpcx/26.3
+    mpirun -np "${n_ranks}" -N "${numproc_node}" --bind-to core --map-by core \
+      "${ARTIFACT}" -f "${DATA}" -pc_type gamg -pc_gamg_square_graph 0 \
+      > "${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=/vol0500/data/ra250029/benchkit-data/petsc-gmres/audikw_1.petscbin
+    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 \
+      > "${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=/home/users/william.dawson/benchkit-data/petsc-gmres/audikw_1.petscbin
+    source /etc/profile.d/modules.sh
+    module load system/ng-dgx nvhpc-hpcx
+    mpirun -np "${n_ranks}" \
+      "${ARTIFACT}" -f "${DATA}" -pc_type gamg -pc_gamg_square_graph 0 \
+      -mat_type aijcusparse \
+      > "${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
+
+solve_time=$(grep "^FOM: ranks=" "${logfile}" | sed -E 's/.*solve_time_s=([0-9.]+).*/\1/')
+
+bk_emit_result \
+  --fom "${solve_time}" \
+  --fom-unit s \
+  --fom-version solve_time \
+  --exp audikw_1 \
+  --nodes "${nodes}" \
+  --numproc-node "${numproc_node}" \
+  --nthreads "${nthreads}" >> "${RESULTS_DIR}/result"
diff --git a/programs/petsc-gmres/src/GMRES-PETSc.cpp b/programs/petsc-gmres/src/GMRES-PETSc.cpp
new file mode 100644
index 0000000..23ff9a4
--- /dev/null
+++ b/programs/petsc-gmres/src/GMRES-PETSc.cpp
@@ -0,0 +1,102 @@
+// 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>
+
+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;
+  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));
+
+  // Figure of merit: wall-clock time of the actual solve (KSPSetUp, which
+  // builds the GAMG hierarchy, runs lazily inside this KSPSolve call, so
+  // setup + iterate are both included -- this is "how long it took to go
+  // from an assembled matrix to a converged solution", independent of
+  // process-launch/link startup cost, which is not part of the solve).
+  PetscLogDouble solve_t0, solve_t1;
+  PetscCall(PetscBarrier((PetscObject)ksp));
+  PetscCall(PetscTime(&solve_t0));
+  PetscCall(KSPSolve(ksp, b, x));
+  PetscCall(PetscBarrier((PetscObject)ksp));
+  PetscCall(PetscTime(&solve_t1));
+  {
+    PetscReal local_dt = (PetscReal)(solve_t1 - solve_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 solve_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;
+}