#!/bin/bash -l
#SBATCH --job-name=smoother-modes
#SBATCH --nodes=36
#SBATCH --ntasks-per-node=8
#SBATCH --cpus-per-task=7
#SBATCH --gpus-per-node=8
#SBATCH --time=1:45:00
#SBATCH --account=phy157_dwf
#SBATCH --gpu-bind=none
#SBATCH --exclusive
#SBATCH --mem=0
#SBATCH -S 0

##############################################################################
# The 1402.2585 p.13 comparison on this machine: adaptive GCR smoothers vs
# the same polynomials frozen (GCRReplaySmoother: recorded for the first
# PolyRecordIters outer steps, then replayed with no inner products), vs a
# Chebyshev 1/x fit on the measured interval.  Same banked point otherwise.
#
#   M1 gcr    / gcr      reference
#   M2 replay / replay   both levels frozen after PolyRecordIters steps
#   M3 replay / gcr      fine frozen only (the fine smoother is the reduction-
#                        heavy one at Nrhs=1)
#   M4 cheb   / gcr      fine Chebyshev [FineChebLo,FineChebHi] order Fso
#   M5 cheb   / cheb
#   M6 gcr    / replay   coarse frozen only
#
# Laptop 8^4 findings (hot config, Ls=4, NBASIS=8): replay/replay converges
# (28 vs 23 outer); cheb on the FINE level diverges there while cheb on the
# coarse level is fine -- the tiny operator has modes the fixed polynomial
# amplifies (left of / off the axis) that GCR handles adaptively.  Production
# spectrum is near-normal with edge 130.5 (shift 0.1) so M4/M5 may behave
# differently; they are cheap to include and cheap to discard.
#
# Readouts per mode: Fouter count, s/RHS, and the FINAL exact-halo residual.
##############################################################################

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=$HOME/ParallelIO/systems/Frontier
source $root/sourceme-rocm7.2.sh

export OMP_NUM_THREADS=7
export FI_MR_CACHE_MONITOR=kdreg2       # REQUIRED for device-buffer MPI on Slingshot: libfabric memhooks monitor (default) is defective, see systems/WorkArounds.txt (libfabric #11451)
export MPICH_GPU_SUPPORT_ENABLED=1
export MPICH_SMP_SINGLE_COPY_MODE=CMA
export MPICH_OFI_NIC_POLICY=GPU
module load libfabric
unset DENSE_GATHER DENSE_GATHER_FORCE DENSE_GATHER_DEBUG DENSE_GATHER_MIN_BYTES

OPTS1="--accelerator-threads 8 --shm 4096 --shm-mpi 1 --device-mem 32000"
vol=48.48.48.96
MPI_GEOM=3.6.4.4

# banked point (2026-08-26 sweep): Css 2.0, Nstep 2, Fso 6, mmax 4, svm 8
export FineSmootherShift=0.1
export FineSmootherOrder=6
export FineSmootherMmax=4
export CoarseSmootherShift=2.0
export CoarseSmootherNstep=2
export CoarseSmootherMmax=2
export CoarseSolverTol=0.05
export CoarseSolverOrder=200
export CoarseSolverMmax=8
export OuterTol=1e-8
export OuterMmax=6
export OuterNstep=12
export SUBSPACE_FILE=/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/subspace_nb64.scidac
unset  SLAB_FILE
export DENSE_SCHUR=1
export DENSE_SCHUR2D=1
export MASS=0.00078
export BLOCK=2.2.3.3
export BLOCK2=4.4.2.4
export L3_TOL=3.0e-1
export L3_MAXIT=2
export L3_NSTEP=50
export DENSE_CC=1
export DENSE_APPLY_PROFILE=1
unset  DENSE_CC_CHECK
export DENSE_SPLITK=128
export DENSE_DEVICE_SUM=4          # cartesian P2P ring ALLGATHER: ~8x fewer bytes than the padded allreduce (=2); no collectives, no size cliff
export SUMMA_HANDSHAKE=1           # ring histogram splits each message into 4-byte handshake (partner wait) + transfer
# SCHUR2D_PROBE: leave unset (probe on by default; =0 disables); probe now also runs 2/3/4 MB
export GRID_ALLOC_NCACHE_LARGE=64
export NRHS=4
export PowerIterations=0
export SmootherCoeffLog=0

# frozen-polynomial controls
export PolyRecordIters=8           # outer steps recorded
export PolyRecordStart=8           # ...starting here: the early-step polynomials are unrepresentative (M3)
export PolyRecordSelect=last       # replay ONE recorded call's polynomial (PB: every individual call beats the coefficient mean)
export PolyRefresh=5               # re-record every 5 outer steps: BFM BfmHDCG.C:2243, k%5==1 -> LdopM1MirsPolyRecord, single call, replayed 4 steps
# Inverse ring-rate hypotheses (2026-08-27): OMP threads and MPI thread level both NIL
# (schur2d_env.job); the cause was partial ring participation -> SCHUR2D_LEAF_SPAN.
export PolyVerbose=1               # frozen smoothers print |r_m|/|r_0| per call: separates 'bad polynomial' from 'linear V-cycle stagnates the outer'
export FineChebLo=3.0              # harvested |R|<0.1 edge / PowerIteration edge x1.05
export FineChebHi=137.0
export CoarseChebLo=8.0
export CoarseChebHi=47.0           # shift 2.0: edge 43.3 x1.08

# Reference: the banked ADAPTIVE optimum, Fso6 / sm4 / Css2.0 / Nstep2 / svm8 ->
# 28.57 s (Nrhs=1), ~14.9 s/RHS (Nrhs=4).  The stationary smoother converged at
# the deliberate overshoot (order 12, fine shift 1.0, coarse Nstep 6); the
# ladder below walks back towards the banked point.  A cell wins if it stays
# convergent AND beats 28.57 s.  Each cell ~5 min.
run_cell () {
  name=$1; export FineSmootherOrder=$2; export FineSmootherShift=$3; export CoarseSmootherNstep=$4
  export FineSmootherMode=$5; export CoarseSmootherMode=$6
  echo "----- $name : Fso=$FineSmootherOrder Fss=$FineSmootherShift Csn=$CoarseSmootherNstep fine=$FineSmootherMode coarse=$CoarseSmootherMode -----"
  fname=log.ladder.$name
  srun -N36 -n288 --kill-on-bad-exit=1 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \
       --mpi ${MPI_GEOM} --grid $vol $OPTS1 --comms-overlap > $fname 2>&1
  echo "  exit $?"; sleep 60
  echo "  $(grep -h 'V2 3-level solve Nrhs' $fname | sed 's/.*V2/V2/' | tr '\n' ' ')"
  echo "  $(grep -h 'Fouter MrhsPGCR: Converged' $fname | sed 's/.*Converged/Converged/' | cut -c1-60 | tr '\n' ' ')"
  echo "  replay per-call |r|/|r0| (Nrhs=1 solve): $(awk '/THREE-level solve, Nrhs = 1/{s=1} s && /Fsmoother replay \|r\|/{v=$NF; n++; t+=v; if(v>mx)mx=v} END{if(n) printf "mean %.4f max %.4f over %d calls", t/n, mx, n}' $fname)"
  grep -h "SCHUR fp64 distributed invert took\|GB/s/rank\|BIG LEAVES\|ring histogram\|^Grid : Message : [0-9.]* s :    >=" $fname | sed 's/^Grid : Message : [0-9.]* s : //' | cut -c1-150 | head -12
}

# Ladder 2 (after L0-L3 of the first ladder): the stationary polynomial
# wants a LARGER shift than the adaptive one -- it must suppress the
# high-mode deviation on average, not per call.  Measured (Csn=2):
#   (12,1.0) 48 outer 31.4 s   (8,1.0) 58 outer 30.9 s   (8,0.5) 60 outer 34.3 s
# Reference (adaptive, Fso6/Fss0.1): 56 outer, 28.57 s.  So walk the SHIFT UP
# at orders 8 and 6, trading shift against order.  Nrhs=1 is the target.
# Ladder 2 result (jobs 5353561/5353660/5356490): stationary polynomials floor
# at ~31 s vs adaptive 28.57 s -- CLOSED.  Cells kept for reference:
#          name          Fso  Fss   Csn  fine    coarse
#run_cell  S1_08_20       8    2.0   2    replay  gcr     # 57 outer 30.55 s
#run_cell  S2_08_30       8    3.0   2    replay  gcr     # 60 / 31.29
#run_cell  S3_06_20       6    2.0   2    replay  gcr     # 69 / 32.74
#run_cell  S4_06_10       6    1.0   2    replay  gcr     # 78 / 35.74
#run_cell  S5_10_20       10   2.0   2    replay  gcr     # 53 / 31.13
#run_cell  S6_coarse6     8    2.0   6    replay  replay  # 61 / 42.53 (intermittent GPU page fault at outer step 22 on one run)
#run_cell  S7_coarse4     8    2.0   4    replay  replay  # 60 / 84.12

# SCHUR2D_LEAF_SPAN sweep (2026-08-27).  The SUMMA ring histogram put 93% of
# the inverse's ring time in 3.7 MB single-block panels: the recursion levels
# below span Pr=18 / Pc=16 run their rings on a few process rows/columns while
# the rest wait.  Big leaves gather each sub-block of <= s blocks to one rank,
# invert locally (rocSOLVER getrf/getri), scatter back.  Readout per cell:
#   "SCHUR fp64 distributed invert took"   -- the number that matters (ref 73-148 s)
#   "BIG LEAVES ... inverse (max over ranks)" -- rocSOLVER time at W = s*480 (unmeasured on a GCD)
#   histogram ">= 2.0 MB" row              -- should shrink to ~nothing at s>=9; handshake column = partner wait
# Solver settings are the banked adaptive point (solve time is the same in every cell; only setup changes).
#          name          Fso  Fss   Csn  fine    coarse
# Sweep result (job 2026-08-27): span 1 132.5 s | 4 37.0 s | 9 27.6 s | 18 30.6 s; VERIFY 3.3939e-5 in all.
# Span 9 optimal: above it the SERIAL leaf chain (32 x 0.53 s at W=4320, 16 x 1.84 s at 8640,
# 287 ranks waiting) dominates; below it the partial rings return.
#export SCHUR2D_LEAF_SPAN=1;  run_cell   I1_span01   6  0.1  2  gcr  gcr     # 132.5 s
#export SCHUR2D_LEAF_SPAN=9;  run_cell   I2_span09   6  0.1  2  gcr  gcr     # 27.6 s
#export SCHUR2D_LEAF_SPAN=18; run_cell   I3_span18   6  0.1  2  gcr  gcr     # 30.6 s
#export SCHUR2D_LEAF_SPAN=4;  run_cell   I4_span04   6  0.1  2  gcr  gcr     # 37.0 s
# Leaf inverse: rocSOLVER getri_batched scales n^1.2-1.8 here (overhead bound).  SCHUR2D_LEAF_LU=1
# replaces it with blocked zgetrf_64 + one identity zgetrs_64 (the N=69120 1-rank path).
# Readout: "BIG LEAVES ... inverse (max over ranks)" vs 0.53 s (span 9) / 1.84 s (span 18);
# if the leaf gets ~3x cheaper, span 18's clean rings (10.3 GB/s) may win overall.
export SCHUR2D_LEAF_LU=1
export SCHUR2D_LEAF_SPAN=9;  run_cell   I5_span09_LU   6  0.1  2  gcr  gcr   # vs I2 27.6 s
export SCHUR2D_LEAF_SPAN=18; run_cell   I6_span18_LU   6  0.1  2  gcr  gcr   # vs I3 30.6 s
unset SCHUR2D_LEAF_SPAN SCHUR2D_LEAF_LU

echo "========================================================="
echo "summary"
for f in log.ladder.I*; do echo "$f: $(grep -h "SCHUR fp64 distributed invert took" $f | sed "s/.*invert took//" | cut -c1-40)  $(grep -h "BIG LEAVES" $f | sed "s/.*leaves,//" | cut -c1-90)"; done
echo "reference inverse (span 1): 73-148 s run to run; solve 28.57 s Nrhs=1"
echo "========================================================="