From 45c1a525b9849dc57b5bd3e04c3320bc3d0c7d85 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 29 Sep 2026 13:47:38 -0400 Subject: [PATCH] Useful script and keep the mixed precision job --- scripts/gcr_polynomial.py | 142 ++++++++++++++ systems/Frontier/pvdagm_mixed_precision.job | 206 ++++++++++++++++++++ 2 files changed, 348 insertions(+) create mode 100644 scripts/gcr_polynomial.py create mode 100644 systems/Frontier/pvdagm_mixed_precision.job diff --git a/scripts/gcr_polynomial.py b/scripts/gcr_polynomial.py new file mode 100644 index 000000000..9d5328011 --- /dev/null +++ b/scripts/gcr_polynomial.py @@ -0,0 +1,142 @@ +#!/usr/bin/env python3 +""" +Reconstruct the GCR smoother polynomial from PGCR coefficient logs +(SmootherCoeffLog=1 in Example_pvdagm_3level_DenseCoarseMatrix, or +LogCoefficients(1) on any PrecGeneralisedConjugateResidualNonHermitian). + +GCR with z = r (trivial preconditioner) builds p_k = P_k(A) r0, +r_k = R_k(A) r0: + + R_0 = P_0 = 1 + R_{k+1}(l) = R_k(l) - a_k l P_k(l) + P_{k+1}(l) = R_{k+1}(l) + sum_j b_{k,j} P_{k-j}(l) + +so after m steps x_m = S(A) r0 with S = sum_k a_k P_k and residual +polynomial R_m = 1 - l S(l), R_m(0) = 1. + +Usage + gcr_polynomial.py LOG [--name Fsmoother] [--lmin 0] [--lmax 3] [--n 600] + writes _poly.dat : l |R_mean(l)| Re R_mean Im R_mean |R|min |R|max |S_mean(l)| + where mean = polynomial built from the per-step MEAN coefficients over + all calls, and min/max are the envelope over the individual calls. + Prints the mean coefficients and their relative spread per step. + gcr_polynomial.py --selftest + runs GCR in numpy on a random non-normal matrix, logs a,b the same way, + and checks R_m(A) r0 == the actual residual to rounding. + +gnuplot: + set logscale y; plot 'Fsmoother_poly.dat' u 1:2 w l t '|R_m| mean', \ + '' u 1:5 w l t 'min', '' u 1:6 w l t 'max' +""" +import re, sys, argparse +import numpy as np +from numpy.polynomial import polynomial as Poly # coefficient arrays, low->high + +A_RE = re.compile(r'(\S+)\s+coeff\[(\d+)\]\s+a=\(([^,]+),([^)]+)\)') +B_RE = re.compile(r'(\S+)\s+coeff\[(\d+)\]((?:\s+b\[\d+\]=\([^)]*\))+)') +BJ_RE = re.compile(r'b\[(\d+)\]=\(([^,]+),([^)]+)\)') + +def parse(path, name): + """-> list of calls; each call = list of (a_k, [b_k0, b_k1, ...]) in step order.""" + calls, cur = [], None + for line in open(path, errors='replace'): + m = A_RE.search(line) + if m and m.group(1) == name: + k = int(m.group(2)); a = complex(float(m.group(3)), float(m.group(4))) + if k == 0: + if cur: calls.append(cur) + cur = [] + if cur is not None: + cur.append([a, []]) + continue + m = B_RE.search(line) + if m and m.group(1) == name and cur is not None: + k = int(m.group(2)) + bs = [complex(float(x[1]), float(x[2])) for x in BJ_RE.findall(m.group(3))] + if k < len(cur): cur[k][1] = bs + if cur: calls.append(cur) + return calls + +def polynomials(coeffs): + """coeffs: list of (a_k, [b_kj]) -> (S, R) as low->high complex coefficient arrays.""" + R = [np.array([1.0+0j])]; P = [np.array([1.0+0j])] + S = np.array([0.0+0j]) + for k, (a, bs) in enumerate(coeffs): + S = Poly.polyadd(S, a*P[k]) + Rn = Poly.polysub(R[k], a*Poly.polymulx(P[k])) # R_{k+1} = R_k - a l P_k + Pn = Rn.copy() + for j, b in enumerate(bs): # P_{k+1} = R_{k+1} + sum b_kj P_{k-j} + Pn = Poly.polyadd(Pn, b*P[k-j]) + R.append(Rn); P.append(Pn) + return S, R[-1] + +def mean_coeffs(calls): + m = min(len(c) for c in calls) + out, spread = [], [] + for k in range(m): + a = np.array([c[k][0] for c in calls]) + nb = min(len(c[k][1]) for c in calls) + bs = [np.array([c[k][1][j] for c in calls]) for j in range(nb)] + out.append((a.mean(), [b.mean() for b in bs])) + sa = a.std()/max(abs(a.mean()),1e-300) + sb = [b.std()/max(abs(b.mean()),1e-300) for b in bs] + spread.append((sa, sb)) + return out, spread + +def selftest(n=40, m=6, mmax=6, seed=1): + rng = np.random.default_rng(seed) + A = np.eye(n)*1.5 + 0.4*rng.standard_normal((n,n)) + 0.2j*rng.standard_normal((n,n)) # non-normal + r0 = rng.standard_normal(n) + 1j*rng.standard_normal(n) + r = r0.copy(); x = np.zeros(n, complex) + p = [r.copy()]; q = [A@r]; qq = [np.vdot(q[0],q[0]).real] + log = [] + for k in range(m): + a = np.vdot(q[k], r)/qq[k] # / + x = x + a*p[k]; r = r - a*q[k] + z = r; Az = A@z; pn = z.copy(); qn = Az.copy(); bs = [] + for back in range(min(k+1, mmax-1)): + b = -(np.vdot(q[k-back], Az)/qq[k-back]).real # real part, as Grid does + pn = pn + b*p[k-back]; qn = qn + b*q[k-back]; bs.append(complex(b)) + p.append(pn); q.append(qn); qq.append(np.vdot(qn,qn).real) + log.append((a, bs)) + S, R = polynomials(log) + def apply(poly, v): + out = np.zeros_like(v); Ak = v.copy() + for c in poly: out = out + c*Ak; Ak = A@Ak + return out + err_r = np.linalg.norm(apply(R, r0) - r)/np.linalg.norm(r) + err_x = np.linalg.norm(apply(S, r0) - x)/np.linalg.norm(x) + print(f"selftest: |R_m(A)r0 - r_m|/|r_m| = {err_r:.2e} |S(A)r0 - x_m|/|x_m| = {err_x:.2e}") + return err_r < 1e-10 and err_x < 1e-10 + +if __name__ == "__main__": + ap = argparse.ArgumentParser() + ap.add_argument("log", nargs="?") + ap.add_argument("--name", default="Fsmoother") + ap.add_argument("--lmin", type=float, default=0.0) + ap.add_argument("--lmax", type=float, default=3.0) + ap.add_argument("--n", type=int, default=600) + ap.add_argument("--selftest", action="store_true") + args = ap.parse_args() + if args.selftest: + sys.exit(0 if selftest() else 1) + calls = parse(args.log, args.name) + if not calls: + sys.exit(f"no '{args.name} coeff[...]' lines found in {args.log}") + mean, spread = mean_coeffs(calls) + print(f"{args.name}: {len(calls)} calls, {len(mean)} steps") + for k,(a,bs) in enumerate(mean): + print(f" step {k}: a = ({a.real:+.5f},{a.imag:+.5f}) rel spread {spread[k][0]:.3f} " + + " ".join(f"b[{j}]={b.real:+.5f} (spread {spread[k][1][j]:.3f})" for j,b in enumerate(bs))) + S, R = polynomials(mean) + print(" S coefficients (low->high):", np.array2string(S, precision=5)) + lam = np.linspace(args.lmin, args.lmax, args.n) + Rm = Poly.polyval(lam, R); Sm = Poly.polyval(lam, S) + per = np.array([np.abs(Poly.polyval(lam, polynomials(c)[1])) for c in calls]) + out = f"{args.name}_poly.dat" + with open(out, "w") as f: + f.write("# lambda |R_mean| ReR ImR |R|min |R|max |S_mean|\n") + for i,l in enumerate(lam): + f.write(f"{l:.6f} {abs(Rm[i]):.6e} {Rm[i].real:.6e} {Rm[i].imag:.6e} " + f"{per[:,i].min():.6e} {per[:,i].max():.6e} {abs(Sm[i]):.6e}\n") + print(f" wrote {out}") diff --git a/systems/Frontier/pvdagm_mixed_precision.job b/systems/Frontier/pvdagm_mixed_precision.job new file mode 100644 index 000000000..c9982cdf2 --- /dev/null +++ b/systems/Frontier/pvdagm_mixed_precision.job @@ -0,0 +1,206 @@ +#!/bin/bash -l +#SBATCH --job-name=pvdagm-mixed-precision +#SBATCH --nodes=36 +#SBATCH --ntasks-per-node=8 +#SBATCH --cpus-per-task=7 +#SBATCH --gpus-per-node=8 +#SBATCH --time=2:00:00 +#SBATCH --account=phy157_dwf +#SBATCH --gpu-bind=none +#SBATCH --exclusive +#SBATCH --mem=0 +#SBATCH -S 0 + +############################################################################## +# The three precision axes of the PVdagM multigrid, on the banked point of +# pvdagm_multigrid.job (48^3x96, Ls=24, nbasis=60, 288 GCDs). The outer +# Krylov is fp64 and exact in every cell, so fp32 below it costs convergence +# rate only, never correctness -- the FINAL exact-halo true residual is the +# proof and is printed by every cell. +# +# 1. fine level inside the preconditioner run time +# 2. coarse + coarse-coarse sector compile -DCOARSE_SINGLE +# 3. distributed dense coarse-coarse invert compile -DGRID_DENSE_INVERSE_SINGLE +# +# The compile-time axes are separate binaries from the same source, all +# plain `make` targets (examples/Makefile.am): +# Example_pvdagm_multigrid fp64 coarse, fp64 dense +# Example_pvdagm_multigrid_fp32coarse fp32 coarse, fp64 dense +# Example_pvdagm_multigrid_fp32coarse_fp32dense fp32 coarse, fp32 dense +# The default basis is 60. The subspace load reads the first 60 vectors +# of the file. +############################################################################## + +LUSTRE=/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle +RUNDIR=$LUSTRE/runs/${SLURM_JOB_NAME}_${SLURM_JOB_ID} +mkdir -p $RUNDIR && cd $RUNDIR && echo "RUNDIR $RUNDIR" + +cat << EOF > select_gpu +#!/bin/bash +export GPU_MAP=(0 1 2 3 7 6 5 4) +export NUMA_MAP=(3 3 1 1 2 2 0 0) +export GPU=\${GPU_MAP[\$SLURM_LOCALID]} +export NUMA=\${NUMA_MAP[\$SLURM_LOCALID]} +export HIP_VISIBLE_DEVICES=\$GPU +unset ROCR_VISIBLE_DEVICES +if [ \$SLURM_PROCID = "0" ]; then echo \$*; fi +exec numactl -m \$NUMA -N \$NUMA \$* +EOF +chmod +x ./select_gpu + +root=/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/MGrewrite/Grid/systems/Frontier +source $root/sourceme-rocm7.2.sh +# sourceme sets FI_HMEM_ROCR_USE_DMABUF=0 -- the standing avoidance of the CXI +# NO_TRANSLATION fault at NRHS>=12 (libfabric #12775). + +export OMP_NUM_THREADS=7 +ulimit -c 0 # no 22 GB GPU core dumps +export FI_MR_CACHE_MONITOR=kdreg2 # site default; device-buffer MPI on Slingshot +export MPICH_GPU_SUPPORT_ENABLED=1 +export MPICH_SMP_SINGLE_COPY_MODE=CMA +export MPICH_OFI_NIC_POLICY=GPU +export GRID_ALLOC_NCACHE_LARGE=64 # allocator large-ring depth, not a solver knob +module load libfabric + +C64D64=$root/examples/Example_pvdagm_multigrid +C32D64=$root/examples/Example_pvdagm_multigrid_fp32coarse +C32D32=$root/examples/Example_pvdagm_multigrid_fp32coarse_fp32dense +# 40000, not 32000: at Nrhs 12 the working set is about 37 GB (outer restart +# history 24.5, fine smoother history 8.2, sources and solutions 4.1), so a +# 32 GB cache evicts. Measured: 808 ms of device-to-host copies in a 20 s +# window, and 9.706 s per rhs against 8.748 with the larger cache. The +# non-evictable total is 6.3 GiB in the fp32 sector, so 40 GB of cache plus +# comms and shm still fits a 64 GB device. +OPTS="--accelerator-threads 8 --shm 4096 --shm-mpi 1 --device-mem 40000 --comms-overlap" +vol=48.48.48.96 +MPI_GEOM=3.6.4.4 + +SUBSPACE=$LUSTRE/subspace_nb64.scidac +CONFIG=/ccs/home/poare/ckpoint_lat.1000 + +############################################################################## +# write_params +# +# Everything else is the banked adaptive operating point of +# pvdagm_multigrid.job, unchanged, so the cells differ only in precision. +############################################################################## +write_params () { + local f=$1 fprec=$2 sloppy=$3 +cat << EOF > $f + + + + 24 + 0.00078 + 1.8 + 1.5 + 0.5 + $CONFIG + 12 + 1 + + + 2233 + 4424 + 9 + $SUBSPACE + $sloppy + $fprec + 0 + + + 0.164 + + + 2.022 + + + 0.052008 + + + 1e-81000612 + + + 9 + + + + +EOF +} + +############################################################################## +# run_cell +############################################################################## +run_cell () { + name=$1; bin=$2; fprec=$3; sloppy=$4 + if [ ! -x "$bin" ]; then + echo "----- $name : SKIPPED, no binary $bin" + echo " build it with a plain make in examples/" + return + fi + write_params params.$name.xml $fprec $sloppy + echo "----- $name : fine=$fprec sloppy=$sloppy bin=$(basename $bin) -----" + export GRID_STDOUT_ROOT=$RUNDIR/mg_$name + srun -N36 -n288 --kill-on-bad-exit=1 ./select_gpu $bin --mpi ${MPI_GEOM} --grid $vol $OPTS \ + --pvdagm-params params.$name.xml --debug-stdout --log Error,Warning,Message,Performance \ + > log.mg.$name 2>&1 + echo " exit $?"; sleep 30 + # Device OOM is reported on stderr as "hipMalloc failed for out of + # memory" and nowhere else, so it must be in this pattern or a cell that + # died of it looks like a cell that printed nothing. + f=$(grep -l "Memory access fault\|NO_TRANSLATION\|hipMalloc failed\|out of memory\|illegal memory access" $GRID_STDOUT_ROOT/*/Grid.stderr.* 2>/dev/null | head -1) + if [ -n "$f" ]; then + echo " FAULT in $f (rank $(basename $f | sed 's/Grid.stderr.//'))"; tail -20 $f | cut -c1-160 + else + echo " certificates, convergence, timing:" + grep -h "Coarse sector precision\|GALERKIN CERTIFICATE\|IMPORT CERTIFICATE\|inversion-source import\|SCHUR .* distributed invert\|VERIFY\|PVdagMMultiGridSolver: Nrhs\|MrhsPGCR: Converged\|FINAL Nrhs .*worst" \ + $GRID_STDOUT_ROOT/0/Grid.stdout.0 2>/dev/null | cut -c1-150 + fi +} + +############################################################################## +# Cells. P0 is the regression gate: with fp64 everywhere it must reproduce +# the banked outer count and residual of pvdagm_multigrid.job. If it does +# not, nothing downstream means anything, and the Galerkin certificate +# (2.3e-15 expected in fp64) separates a coarsening change from a solve one. +############################################################################## + +# name binary fineprec sloppy +run_cell P0_fp64 $C64D64 fp64 1 +run_cell P1_fp32fine $C64D64 fp32 1 +run_cell P2_fp32coarse $C32D64 fp64 1 +run_cell P3_fp32both $C32D64 fp32 1 +run_cell P4_denseF32 $C32D32 fp32 1 + +# Halo wire format is orthogonal to arithmetic: on an fp32 operator +# FineSloppyComms is bf16 compression, on an fp64 operator it is fp32. One +# exact-halo cell to price it inside the preconditioner. +run_cell P5_exacthalo $C32D64 fp32 0 + +############################################################################## +# What to read. +# +# P0: the gate. Outer count and s/RHS against the banked numbers. +# +# P1 vs P0: fp32 fine level. Laptop 8^4 saw the same outer count and ~12% +# less wall; here the fine level is the dominant cost and the halo is real, +# so this is where the fp32 case is actually made or lost. +# +# P2 vs P0: fp32 coarse sector. Watch the Galerkin certificate move from +# ~2e-15 to ~3e-7 -- that is fp32 rounding in the coarsening, not an error -- +# and check the outer count does not move with it. Halves the coarse +# operator's storage and its comms. +# +# P4 vs P3: fp32 dense coarse-coarse inversion. The reading is the VERIFY +# certificate ||A Ainv x - x||/||x||, which scales as kappa(A_cc) * eps32 and +# so measures the conditioning of the coarse-coarse operator at production +# size. The laptop's N=128 was too small to say anything (5.3e-6 -> 6.7e-6); +# at N=69120 it is a real test. If VERIFY degrades but the outer count holds, +# the fp32 factorisation is fine and halves that footprint; if the outer count +# moves, the fallback is fp32 factorise plus one fp64 refinement step, which +# is two extra dense applies and is not yet written. +# +# P5 vs P3: what the reduced-precision halo is worth once the arithmetic is +# already fp32 (bf16 wire on an fp32 operator). +##############################################################################