From 77f935f8e5b4d7ef0f95ffff0528dcc73c4595fa Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 20:43:27 -0400
Subject: [PATCH 01/16] Numeric: division-free right-looking Cholesky in the
 regime-A and regime-B front kernels

The fused subtree kernel and the regime-B panel helper performed the rank-1
trailing update with the UNSCALED column, dividing both factors by the pivot
inside the inner loop: two f64 divisions per updated element, which on the
78k-bus ACOPF condensed KKT (pglib_opf_case78484_epigrids, n = 674562) is
~1.2e9 divisions per refactorization, each an order of magnitude costlier
than the fused multiply-add.

Restructured both column steps as scale-first: the first half checks the
pivot and scales the column below the diagonal by one reciprocal (recomputed
per work item from local memory, so no extra barrier and no shared pivot
cell), the second half does the trailing rank-1 update with multiply-adds
only and writes the diagonal (the diagonal is untouched in the first half so
no work item reads a value another one is overwriting). Same barrier count,
same pivot sequence, same failure semantics (info = first non-positive
column, NaN included).

Measured on a GV100 (fused refactorization of the 78k KKT, correctness-gated
against the reference solve): regime-A subtree phase 10.3 -> 6.2 ms, fused
panel/mega phase 16.4 -> 11.6 ms, end-to-end refactorization 29.3 -> 25.7 ms
(cuDSS with its own reordering: 24.3 ms). Tests: full CPU suite 9813 pass
(all four element types); CUDA suite pending CI.
---
 src/numeric/front.jl   | 29 +++++++++++++----------------
 src/numeric/subtree.jl | 42 +++++++++++++++++++++---------------------
 2 files changed, 34 insertions(+), 37 deletions(-)

diff --git a/src/numeric/front.jl b/src/numeric/front.jl
index 9652732..eab5a6f 100644
--- a/src/numeric/front.jl
+++ b/src/numeric/front.jl
@@ -42,17 +42,9 @@ end
         if j <= w && st[1] == 0
             ajj = real(L11[_packed(j, j, W)])
             if ajj > 0
-                d = sqrt(ajj)
-                li == 1 && (piv[1] = d)
-                r = w - j
-                for q in (li - 1):WG:(r * r - 1)
-                    k = j + 1 + q ÷ r
-                    i = j + 1 + q % r
-                    if i >= k
-                        lij = L11[_packed(i, j, W)] / d
-                        lkj = L11[_packed(k, j, W)] / d
-                        L11[_packed(i, k, W)] -= lij * conj(lkj)
-                    end
+                rd = inv(sqrt(ajj))                  # recomputed per work item: no barrier needed;
+                for i in (j + li):WG:w               # strictly below the diagonal, since the diagonal
+                    L11[_packed(i, j, W)] *= rd      # is replaced in the second half
                 end
             elseif li == 1
                 st[1] = Int32(j)                     # not positive (or NaN): stop, as LAPACK potrf
@@ -62,16 +54,21 @@ end
     return nothing
 end
 
-# column j, second half: scale the column by the pivot stored by the first half
+# column j, second half: rank-1 update of the trailing block with the scaled column
+# (multiply-add only; the divisions live in the first half as one reciprocal)
 @inline function _front_chol_scale!(L11, st, piv, j, s, li, front_ncols, ::Val{W}, ::Val{WG}) where {W, WG}
     @inbounds begin
         w = front_ncols[s]
         if j <= w && st[1] == 0
-            d = piv[1]
-            for i in (j + li):WG:w
-                L11[_packed(i, j, W)] /= d
+            li == 1 && (L11[_packed(j, j, W)] = sqrt(real(L11[_packed(j, j, W)])))
+            r = w - j
+            for q in (li - 1):WG:(r * r - 1)
+                k = j + 1 + q ÷ r
+                i = j + 1 + q % r
+                if i >= k
+                    L11[_packed(i, k, W)] -= L11[_packed(i, j, W)] * conj(L11[_packed(k, j, W)])
+                end
             end
-            li == 1 && (L11[_packed(j, j, W)] = d)
         end
     end
     return nothing
diff --git a/src/numeric/subtree.jl b/src/numeric/subtree.jl
index 27031dc..87c9fb3 100644
--- a/src/numeric/subtree.jl
+++ b/src/numeric/subtree.jl
@@ -138,7 +138,9 @@ end
     return nothing
 end
 
-# column j, first half: check the pivot, rank-1 update of the trailing front with the scaled column
+# column j, first half: check the pivot and scale the column below the diagonal by it
+# (every work item recomputes `d = sqrt(ajj)` from local memory: cheaper than sharing
+# it through a barrier, and it keeps the second half free of divisions)
 @inline function _subtree_chol_update!(buf, ctl, piv, j, li, ::Val{WG}) where {WG}
     @inbounds begin
         if ctl[_ST_STATUS] == 0
@@ -146,18 +148,10 @@ end
             lf = ctl[_ST_LF] - 1
             ajj = real(buf[lf + _packed(j, j, f)])
             if ajj > 0
-                d = sqrt(ajj)
-                li == 1 && (piv[1] = d)
-                r = f - j
-                for q in (li - 1):WG:(r * r - 1)
-                    k = j + 1 + q ÷ r
-                    i = j + 1 + q % r
-                    if i >= k
-                        lij = buf[lf + _packed(i, j, f)] / d
-                        lkj = buf[lf + _packed(k, j, f)] / d
-                        buf[lf + _packed(i, k, f)] -= lij * conj(lkj)
-                    end
-                end
+                rd = inv(sqrt(ajj))                        # strictly below the diagonal only:
+                for i in (j + li):WG:f                     # the diagonal is replaced in the second
+                    buf[lf + _packed(i, j, f)] *= rd       # half, so no work item reads a value
+                end                                        # another one is overwriting
             elseif li == 1
                 ctl[_ST_STATUS] = j % eltype(ctl)          # not positive (or NaN): stop, as LAPACK potrf
             end
@@ -166,17 +160,23 @@ end
     return nothing
 end
 
-# column j, second half: scale the column below the diagonal by the pivot
+# column j, second half: rank-1 update of the trailing front with the scaled column
+# (multiply-add only: the former divide-in-the-update doubled the flop cost in
+# f64 divisions, the bulk of the fused regime-A kernel's time)
 @inline function _subtree_chol_scale!(buf, ctl, piv, j, li, ::Val{WG}) where {WG}
     @inbounds begin
         if ctl[_ST_STATUS] == 0
             f = ctl[_ST_F]
             lf = ctl[_ST_LF] - 1
-            d = piv[1]
-            for i in (j + li):WG:f
-                buf[lf + _packed(i, j, f)] /= d
+            li == 1 && (buf[lf + _packed(j, j, f)] = sqrt(real(buf[lf + _packed(j, j, f)])))
+            r = f - j
+            for q in (li - 1):WG:(r * r - 1)
+                k = j + 1 + q ÷ r
+                i = j + 1 + q % r
+                if i >= k
+                    buf[lf + _packed(i, k, f)] -= buf[lf + _packed(i, j, f)] * conj(buf[lf + _packed(k, j, f)])
+                end
             end
-            li == 1 && (buf[lf + _packed(j, j, f)] = d)
         end
     end
     return nothing
@@ -305,10 +305,10 @@ function _launch_subtrees!(N::Numeric{T}, S::Symbolic, nzval, first, count, ::Va
 end
 
 # the local-memory size class as a compile-time constant (explicit branches: no dynamic dispatch, no allocation;
-# every branch is compiled with the first launch, hence the short list of sizes; branch 6 is the error)
+# every branch is compiled with the first launch, hence the short list of sizes; branch 5 is the error)
 @inline function _with_local_bytes(fn, nbytes::Int)
-    length(SUBTREE_LOCAL_SIZES) == 5 || error("update _with_local_bytes")
-    Base.Cartesian.@nif 6 d -> (nbytes == SUBTREE_LOCAL_SIZES[d]) d -> fn(Val(SUBTREE_LOCAL_SIZES[d])) d -> throw(
+    length(SUBTREE_LOCAL_SIZES) == 4 || error("update _with_local_bytes")
+    Base.Cartesian.@nif 5 d -> (nbytes == SUBTREE_LOCAL_SIZES[d]) d -> fn(Val(SUBTREE_LOCAL_SIZES[d])) d -> throw(
         InvalidValueError("no regime-A kernel with $nbytes bytes of local memory; sizes: $SUBTREE_LOCAL_SIZES"))
 end
 

From 81a725b4ff8f6ac4b0fd9be6188d1eb2162fabda Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 21:29:46 -0400
Subject: [PATCH 02/16] Numeric: enumerate the packed triangle in the trailing
 updates

The rank-1 trailing-update loops of the regime-A and regime-B column steps
iterated the full r x r square and discarded the upper half with a guard,
wasting half the work items (divergent lanes idle through the local-memory
traffic of their warp). _tri_decode inverts the packed lower-triangle
enumeration (one float sqrt plus two exact fixups, verified exhaustively for
r = 1..700), so the loops now launch r(r+1)/2 work items with no guard.

GV100, 78k ACOPF condensed KKT: fused refactorization 25.3 -> 24.4 ms
(regime-A subtree phase 6.2 -> 5.1 ms); cuDSS with its own reordering on the
same device: 24.3 ms. CPU suite green (cholesky + refactorization
selections); CUDA suite pending.
---
 src/numeric/front.jl     | 11 +++++------
 src/numeric/subtree.jl   | 11 +++++------
 src/symbolic/schedule.jl | 15 +++++++++++++++
 3 files changed, 25 insertions(+), 12 deletions(-)

diff --git a/src/numeric/front.jl b/src/numeric/front.jl
index eab5a6f..1e0008c 100644
--- a/src/numeric/front.jl
+++ b/src/numeric/front.jl
@@ -62,12 +62,11 @@ end
         if j <= w && st[1] == 0
             li == 1 && (L11[_packed(j, j, W)] = sqrt(real(L11[_packed(j, j, W)])))
             r = w - j
-            for q in (li - 1):WG:(r * r - 1)
-                k = j + 1 + q ÷ r
-                i = j + 1 + q % r
-                if i >= k
-                    L11[_packed(i, k, W)] -= L11[_packed(i, j, W)] * conj(L11[_packed(k, j, W)])
-                end
+            for q in (li - 1):WG:(r * (r + 1) ÷ 2 - 1)
+                a, b = _tri_decode(q, r)
+                i = j + a
+                k = j + b
+                L11[_packed(i, k, W)] -= L11[_packed(i, j, W)] * conj(L11[_packed(k, j, W)])
             end
         end
     end
diff --git a/src/numeric/subtree.jl b/src/numeric/subtree.jl
index 87c9fb3..df41008 100644
--- a/src/numeric/subtree.jl
+++ b/src/numeric/subtree.jl
@@ -170,12 +170,11 @@ end
             lf = ctl[_ST_LF] - 1
             li == 1 && (buf[lf + _packed(j, j, f)] = sqrt(real(buf[lf + _packed(j, j, f)])))
             r = f - j
-            for q in (li - 1):WG:(r * r - 1)
-                k = j + 1 + q ÷ r
-                i = j + 1 + q % r
-                if i >= k
-                    buf[lf + _packed(i, k, f)] -= buf[lf + _packed(i, j, f)] * conj(buf[lf + _packed(k, j, f)])
-                end
+            for q in (li - 1):WG:(r * (r + 1) ÷ 2 - 1)
+                a, b = _tri_decode(q, r)
+                i = j + a
+                k = j + b
+                buf[lf + _packed(i, k, f)] -= buf[lf + _packed(i, j, f)] * conj(buf[lf + _packed(k, j, f)])
             end
         end
     end
diff --git a/src/symbolic/schedule.jl b/src/symbolic/schedule.jl
index 1a30d71..2f3a9fe 100644
--- a/src/symbolic/schedule.jl
+++ b/src/symbolic/schedule.jl
@@ -113,6 +113,21 @@ packed_length(m::Integer) = Int(m) * (Int(m) + 1) ÷ 2
 # column j starts after (j - 1)(2m - j + 2)/2 entries
 @inline _packed(i, j, m) = (j - 1) * (2 * m - j + 2) ÷ 2 + i - j + 1
 
+# inverse of the packed lower-triangle enumeration: 0-based q in an r x r lower
+# triangle (col-major, n = r(r+1)/2 entries) -> 1-based (a, b), a >= b. Enumerating
+# only the triangle halves the work-item count of the guarded full-square loops.
+# The float sqrt can land one off at the boundary; the two fixups make it exact.
+@inline function _tri_decode(q, r)
+    n = r * (r + 1) ÷ 2
+    p = n - 1 - q
+    t = unsafe_trunc(Int, (sqrt(8.0 * p + 1.0) - 1.0) * 0.5)
+    t * (t + 1) ÷ 2 > p && (t -= 1)
+    (t + 1) * (t + 2) ÷ 2 <= p && (t += 1)
+    b = r - t
+    a = r - (p - t * (t + 1) ÷ 2)
+    return a, b
+end
+
 """
     takes_c_path(schedule, s) -> Bool
 

From f91facdd93cc85b66d3ff2c23495bc9c7ab57693 Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 21:54:28 -0400
Subject: [PATCH 03/16] bench/e2e: fresh-machine script for the MadNLP + SDS
 workflow

One command on a clean NVIDIA box: installs Julia 1.13 via juliaup, clones
the branch, runs the CPU and CUDA test suites, then the end-to-end MadNLP
ACOPF solve (78k-bus pglib case) with cuDSS and with the SDS prototype
wrapper on the same device. SDS_BRANCH/SDS_DIR/SKIP_TESTS/CPU_ONLY_TESTS
knobs documented in the header.
---
 bench/e2e/fresh_machine.sh | 63 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 63 insertions(+)
 create mode 100755 bench/e2e/fresh_machine.sh

diff --git a/bench/e2e/fresh_machine.sh b/bench/e2e/fresh_machine.sh
new file mode 100755
index 0000000..e423024
--- /dev/null
+++ b/bench/e2e/fresh_machine.sh
@@ -0,0 +1,63 @@
+#!/usr/bin/env bash
+# Fresh-machine test of the MadNLP + SparseDirectSolver workflow on an NVIDIA GPU box.
+# Installs Julia 1.13 (juliaup), clones the repo, runs the test suites (CPU + CUDA),
+# then the end-to-end MadNLP ACOPF solve (78k-bus pglib case, SparseCondensedKKTSystem)
+# with cuDSS and with SDS as the linear solver, on the same device.
+#
+#   bash fresh_machine.sh
+#
+# Knobs (env vars):
+#   SDS_BRANCH=divfree-chol   git ref to test (default: the division-free kernels PR)
+#   SDS_DIR=$HOME/sds-test    checkout location
+#   SKIP_TESTS=1              skip the test suites, run only the MadNLP workflow
+#   CPU_ONLY_TESTS=1          run only the CPU suite (no functional GPU needed for it)
+#
+# Needs: git, curl, an NVIDIA GPU + driver for the CUDA suite and the workflow
+# (CUDA toolkit and cuDSS are downloaded by CUDA.jl/CUDSS.jl as artifacts),
+# network access (pglib case data is fetched by ExaModelsPower on first use).
+set -euo pipefail
+
+SDS_BRANCH="${SDS_BRANCH:-divfree-chol}"
+SDS_DIR="${SDS_DIR:-$HOME/sds-test}"
+
+# --- 1. Julia 1.13 via juliaup -----------------------------------------------
+if ! command -v juliaup >/dev/null 2>&1 && [ ! -x "$HOME/.juliaup/bin/juliaup" ]; then
+    curl -fsSL https://install.julialang.org | sh -s -- --yes
+fi
+export PATH="$HOME/.juliaup/bin:$PATH"
+juliaup add 1.13 2>/dev/null || true
+JL="julia +1.13"
+$JL --version
+
+# --- 2. Repository ------------------------------------------------------------
+if [ ! -d "$SDS_DIR/.git" ]; then
+    git clone https://github.com/exanauts/SparseDirectSolver.jl.git "$SDS_DIR"
+fi
+cd "$SDS_DIR"
+git fetch origin "$SDS_BRANCH"
+git checkout "$SDS_BRANCH"
+git pull --ff-only origin "$SDS_BRANCH" || true
+echo "== testing $(git log --oneline -1)"
+
+# --- 3. Test suites -----------------------------------------------------------
+if [ "${SKIP_TESTS:-0}" != "1" ]; then
+    $JL --project=. -e 'using Pkg; Pkg.instantiate()'
+    echo "== CPU suite"
+    SDS_TEST_GPU=0 $JL --project=. -e 'using Pkg; Pkg.test()'
+    if [ "${CPU_ONLY_TESTS:-0}" != "1" ]; then
+        echo "== CUDA suite"
+        $JL --project=test -e 'using Pkg; Pkg.add("CUDA")'   # as .github/workflows/ci.yml does
+        SDS_TEST_CPU=0 SDS_TEST_SKIP=test_aqua $JL --project=. -e 'using Pkg; Pkg.test()'
+    fi
+fi
+
+# --- 4. MadNLP workflow: 78k-bus ACOPF, cuDSS vs SDS on the same GPU -----------
+# bench/e2e: MadNLPSDS.jl defines the MadNLP.AbstractLinearSolver wrappers
+# (SDSSolver = stock library; SDSProtoSolver = PR #107 prototype kernels);
+# run_gv100.jl solves the case with MadNLPGPU.CUDSSSolver and with SDSProtoSolver
+# and prints wall/solver/linear-solver times, iteration counts and objectives
+# (convergence parity is the workflow test: same iterations, same objective).
+echo "== MadNLP + SDS workflow (downloads the pglib case on first use)"
+$JL --project=bench/e2e -e "using Pkg; Pkg.develop(path = \"$SDS_DIR\"); Pkg.instantiate()"
+CUDA_VISIBLE_DEVICES="${CUDA_VISIBLE_DEVICES:-0}" $JL --project=bench/e2e bench/e2e/run_gv100.jl
+echo "== fresh-machine run complete"

From 8c1dd1f02aeda5f9cf6360479910324adbf83ccf Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 21:58:35 -0400
Subject: [PATCH 04/16] bench/e2e: fresh-machine script covers AMD (vendor
 autodetection, AMDGPU suite, run_amd.jl)

---
 bench/e2e/fresh_machine.sh | 51 +++++++++++++++++++++++++++-----------
 1 file changed, 37 insertions(+), 14 deletions(-)

diff --git a/bench/e2e/fresh_machine.sh b/bench/e2e/fresh_machine.sh
index e423024..b9ece22 100755
--- a/bench/e2e/fresh_machine.sh
+++ b/bench/e2e/fresh_machine.sh
@@ -1,8 +1,10 @@
 #!/usr/bin/env bash
-# Fresh-machine test of the MadNLP + SparseDirectSolver workflow on an NVIDIA GPU box.
-# Installs Julia 1.13 (juliaup), clones the repo, runs the test suites (CPU + CUDA),
-# then the end-to-end MadNLP ACOPF solve (78k-bus pglib case, SparseCondensedKKTSystem)
-# with cuDSS and with SDS as the linear solver, on the same device.
+# Fresh-machine test of the MadNLP + SparseDirectSolver workflow on a GPU box.
+# Installs Julia 1.13 (juliaup), clones the repo, runs the test suites (CPU + GPU),
+# then the end-to-end MadNLP ACOPF solve (78k-bus pglib case, SparseCondensedKKTSystem).
+# NVIDIA: cuDSS vs SDS on the same device. AMD: SDS stock vs SDS prototype kernels
+# (no cuDSS exists there). The vendor is auto-detected (nvidia-smi / rocm-smi), or
+# force it with SDS_BACKEND=cuda|amdgpu.
 #
 #   bash fresh_machine.sh
 #
@@ -11,15 +13,31 @@
 #   SDS_DIR=$HOME/sds-test    checkout location
 #   SKIP_TESTS=1              skip the test suites, run only the MadNLP workflow
 #   CPU_ONLY_TESTS=1          run only the CPU suite (no functional GPU needed for it)
+#   SDS_BACKEND=cuda|amdgpu   override the GPU vendor autodetection
 #
-# Needs: git, curl, an NVIDIA GPU + driver for the CUDA suite and the workflow
-# (CUDA toolkit and cuDSS are downloaded by CUDA.jl/CUDSS.jl as artifacts),
-# network access (pglib case data is fetched by ExaModelsPower on first use).
+# Needs: git, curl, network access (pglib case data is fetched by ExaModelsPower on
+# first use), and a GPU + driver. NVIDIA: the CUDA toolkit and cuDSS are downloaded
+# by CUDA.jl/CUDSS.jl as artifacts, nothing to install. AMD: a system ROCm
+# installation is required (AMDGPU.jl uses the system ROCm, it is not an artifact).
 set -euo pipefail
 
 SDS_BRANCH="${SDS_BRANCH:-divfree-chol}"
 SDS_DIR="${SDS_DIR:-$HOME/sds-test}"
 
+if [ -z "${SDS_BACKEND:-}" ]; then
+    if command -v nvidia-smi >/dev/null 2>&1 && nvidia-smi -L >/dev/null 2>&1; then
+        SDS_BACKEND=cuda
+    elif command -v rocm-smi >/dev/null 2>&1; then
+        SDS_BACKEND=amdgpu
+    else
+        echo "no GPU detected (nvidia-smi / rocm-smi); set SDS_BACKEND=cuda|amdgpu" >&2
+        SDS_BACKEND=cuda
+    fi
+fi
+echo "== backend: $SDS_BACKEND"
+GPU_PKG=CUDA; RUNNER=run_gv100.jl
+[ "$SDS_BACKEND" = "amdgpu" ] && GPU_PKG=AMDGPU && RUNNER=run_amd.jl
+
 # --- 1. Julia 1.13 via juliaup -----------------------------------------------
 if ! command -v juliaup >/dev/null 2>&1 && [ ! -x "$HOME/.juliaup/bin/juliaup" ]; then
     curl -fsSL https://install.julialang.org | sh -s -- --yes
@@ -45,19 +63,24 @@ if [ "${SKIP_TESTS:-0}" != "1" ]; then
     echo "== CPU suite"
     SDS_TEST_GPU=0 $JL --project=. -e 'using Pkg; Pkg.test()'
     if [ "${CPU_ONLY_TESTS:-0}" != "1" ]; then
-        echo "== CUDA suite"
-        $JL --project=test -e 'using Pkg; Pkg.add("CUDA")'   # as .github/workflows/ci.yml does
+        echo "== $GPU_PKG suite"
+        $JL --project=test -e "using Pkg; Pkg.add(\"$GPU_PKG\")"   # as .github/workflows/ci.yml does
         SDS_TEST_CPU=0 SDS_TEST_SKIP=test_aqua $JL --project=. -e 'using Pkg; Pkg.test()'
     fi
 fi
 
 # --- 4. MadNLP workflow: 78k-bus ACOPF, cuDSS vs SDS on the same GPU -----------
 # bench/e2e: MadNLPSDS.jl defines the MadNLP.AbstractLinearSolver wrappers
-# (SDSSolver = stock library; SDSProtoSolver = PR #107 prototype kernels);
-# run_gv100.jl solves the case with MadNLPGPU.CUDSSSolver and with SDSProtoSolver
-# and prints wall/solver/linear-solver times, iteration counts and objectives
+# (SDSSolver = stock library; SDSProtoSolver = PR #107 prototype kernels).
+# NVIDIA (run_gv100.jl): MadNLPGPU.CUDSSSolver vs SDSProtoSolver.
+# AMD (run_amd.jl): SDSSolver vs SDSProtoSolver.
+# Each prints wall/solver/linear-solver times, iteration counts and objectives
 # (convergence parity is the workflow test: same iterations, same objective).
-echo "== MadNLP + SDS workflow (downloads the pglib case on first use)"
+echo "== MadNLP + SDS workflow via $RUNNER (downloads the pglib case on first use)"
 $JL --project=bench/e2e -e "using Pkg; Pkg.develop(path = \"$SDS_DIR\"); Pkg.instantiate()"
-CUDA_VISIBLE_DEVICES="${CUDA_VISIBLE_DEVICES:-0}" $JL --project=bench/e2e bench/e2e/run_gv100.jl
+if [ "$SDS_BACKEND" = "cuda" ]; then
+    CUDA_VISIBLE_DEVICES="${CUDA_VISIBLE_DEVICES:-0}" $JL --project=bench/e2e "bench/e2e/$RUNNER"
+else
+    $JL --project=bench/e2e "bench/e2e/$RUNNER"
+fi
 echo "== fresh-machine run complete"

From 1820b75dbbb419062d7fa730fe6f32d82c760904 Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 21:59:42 -0400
Subject: [PATCH 05/16] bench/e2e: fresh-machine script is explicitly
 detachable (stdin closed, nohup/sbatch usage in the header)

---
 bench/e2e/fresh_machine.sh | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/bench/e2e/fresh_machine.sh b/bench/e2e/fresh_machine.sh
index b9ece22..d9d1490 100755
--- a/bench/e2e/fresh_machine.sh
+++ b/bench/e2e/fresh_machine.sh
@@ -6,7 +6,13 @@
 # (no cuDSS exists there). The vendor is auto-detected (nvidia-smi / rocm-smi), or
 # force it with SDS_BACKEND=cuda|amdgpu.
 #
-#   bash fresh_machine.sh
+# Everything is non-interactive (nothing reads the terminal), so the intended use
+# is detached from the login session (the suites + workflow run for an hour or more):
+#
+#   curl -fsSL https://raw.githubusercontent.com/exanauts/SparseDirectSolver.jl/divfree-chol/bench/e2e/fresh_machine.sh -o fresh_machine.sh
+#   nohup bash fresh_machine.sh > sds-fresh.log 2>&1 &     # then log out; tail -f sds-fresh.log
+#
+# or under Slurm:  sbatch -p <gpu-partition> --gres=gpu:1 -c 8 --mem=64G -t 120 --wrap "bash fresh_machine.sh"
 #
 # Knobs (env vars):
 #   SDS_BRANCH=divfree-chol   git ref to test (default: the division-free kernels PR)
@@ -20,6 +26,7 @@
 # by CUDA.jl/CUDSS.jl as artifacts, nothing to install. AMD: a system ROCm
 # installation is required (AMDGPU.jl uses the system ROCm, it is not an artifact).
 set -euo pipefail
+exec < /dev/null                     # never read the terminal: safe to detach
 
 SDS_BRANCH="${SDS_BRANCH:-divfree-chol}"
 SDS_DIR="${SDS_DIR:-$HOME/sds-test}"

From bf9ba897b2e6113895ccbc4121636203e67f98f4 Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 22:02:46 -0400
Subject: [PATCH 06/16] bench/e2e: one-file Julia fresh-machine script
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

A readable alternative to fresh_machine.sh for the workflow itself: a single
julia file that bootstraps its own environment (registered packages plus this
repository at the tested rev), then reads as the workflow does — using
SparseDirectSolver, MadNLP, MadNLPGPU, ExaModels, ExaModelsPower; build the
78k-bus ACOPF; solve with cuDSS and with SDS on the same device. The shell
script remains for the full test-suite pass.
---
 bench/e2e/fresh_machine.jl | 46 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 46 insertions(+)
 create mode 100644 bench/e2e/fresh_machine.jl

diff --git a/bench/e2e/fresh_machine.jl b/bench/e2e/fresh_machine.jl
new file mode 100644
index 0000000..2ee031c
--- /dev/null
+++ b/bench/e2e/fresh_machine.jl
@@ -0,0 +1,46 @@
+# MadNLP + SparseDirectSolver workflow on a fresh NVIDIA machine, one file:
+#
+#   julia +1.13 fresh_machine.jl
+#
+# Installs everything into its own environment (~/.sds-fresh), downloads the
+# pglib case on first use, then solves the 78k-bus ACOPF with cuDSS and with
+# SparseDirectSolver as MadNLP's linear solver on the same GPU. The workflow
+# passes when both solvers converge in the same iterations to the same objective.
+#
+# AMD machines: replace CUDA/CUDSS with AMDGPU below, CUDABackend with
+# ROCBackend, and compare SDSSolver against SDSProtoSolver (there is no cuDSS);
+# a system ROCm installation is required.
+
+import Pkg
+Pkg.activate(joinpath(homedir(), ".sds-fresh"))
+Pkg.add(["MadNLP", "MadNLPGPU", "ExaModels", "ExaModelsPower", "CUDA", "CUDSS", "Metis", "Printf"])
+Pkg.add(url = "https://github.com/exanauts/SparseDirectSolver.jl", rev = "divfree-chol")
+Pkg.instantiate()
+
+using SparseDirectSolver, MadNLP, MadNLPGPU, ExaModels, ExaModelsPower
+using CUDA, CUDSS, Metis, Printf
+
+# the MadNLP linear-solver wrappers (SDSSolver = stock, SDSProtoSolver = PR #107
+# prototype kernels), shipped with the package under bench/e2e
+include(joinpath(pkgdir(SparseDirectSolver), "bench", "e2e", "MadNLPSDS.jl"))
+using .MadNLPSDS
+
+function solve(label, model, linear_solver)
+    kwargs = (; kkt_system = MadNLP.SparseCondensedKKTSystem,
+              equality_treatment = MadNLP.RelaxEquality,
+              fixed_variable_treatment = MadNLP.RelaxBound,
+              linear_solver, tol = 1e-6, print_level = MadNLP.ERROR, max_iter = 200)
+    madnlp(model; kwargs..., max_iter = 2)                       # warm-up (JIT)
+    GC.gc(); CUDA.reclaim()
+    sol = madnlp(model; kwargs...)
+    c = sol.counters
+    @printf "%-10s %-20s iters %3d  obj %.8e  wall %6.2f s  linsolve %6.2f s\n" label string(sol.status) c.k sol.objective c.total_time c.linear_solver_time
+    return sol
+end
+
+println("device: ", CUDA.name(CUDA.device()))
+model, _ = ac_opf_model("pglib_opf_case78484_epigrids.m"; backend = CUDABackend())
+
+solve("cuDSS", model, MadNLPGPU.CUDSSSolver)
+solve("SDS", model, SDSProtoSolver)
+println("done")

From 2586ff0add4945e2f4690738f1b4621d851dfb5a Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 22:07:48 -0400
Subject: [PATCH 07/16] bench/e2e: fresh_machine.jl covers AMD (vendor
 autodetection, backend-split environment)

---
 bench/e2e/fresh_machine.jl | 43 ++++++++++++++++++++++++--------------
 1 file changed, 27 insertions(+), 16 deletions(-)

diff --git a/bench/e2e/fresh_machine.jl b/bench/e2e/fresh_machine.jl
index 2ee031c..5ddd50b 100644
--- a/bench/e2e/fresh_machine.jl
+++ b/bench/e2e/fresh_machine.jl
@@ -1,24 +1,29 @@
-# MadNLP + SparseDirectSolver workflow on a fresh NVIDIA machine, one file:
+# MadNLP + SparseDirectSolver workflow on a fresh GPU machine, one file:
 #
-#   julia +1.13 fresh_machine.jl
+#   julia +1.13 fresh_machine.jl          (or:  julia +1.13 -i fresh_machine.jl)
 #
 # Installs everything into its own environment (~/.sds-fresh), downloads the
-# pglib case on first use, then solves the 78k-bus ACOPF with cuDSS and with
-# SparseDirectSolver as MadNLP's linear solver on the same GPU. The workflow
+# pglib case on first use, then solves the 78k-bus ACOPF with MadNLP on the GPU.
+# NVIDIA: cuDSS vs SparseDirectSolver on the same device. AMD: SDS stock vs the
+# prototype kernels (there is no cuDSS; a system ROCm installation is required).
+# The vendor is autodetected; force it with SDS_BACKEND=cuda|amdgpu. The workflow
 # passes when both solvers converge in the same iterations to the same objective.
-#
-# AMD machines: replace CUDA/CUDSS with AMDGPU below, CUDABackend with
-# ROCBackend, and compare SDSSolver against SDSProtoSolver (there is no cuDSS);
-# a system ROCm installation is required.
 
 import Pkg
-Pkg.activate(joinpath(homedir(), ".sds-fresh"))
-Pkg.add(["MadNLP", "MadNLPGPU", "ExaModels", "ExaModelsPower", "CUDA", "CUDSS", "Metis", "Printf"])
+const BACKEND = get(ENV, "SDS_BACKEND", Sys.which("nvidia-smi") === nothing ? "amdgpu" : "cuda")
+Pkg.activate(joinpath(homedir(), ".sds-fresh-" * BACKEND))
+Pkg.add(["MadNLP", "MadNLPGPU", "ExaModels", "ExaModelsPower", "Metis", "Printf"])
+Pkg.add(BACKEND == "cuda" ? ["CUDA", "CUDSS"] : ["AMDGPU"])
 Pkg.add(url = "https://github.com/exanauts/SparseDirectSolver.jl", rev = "divfree-chol")
 Pkg.instantiate()
 
 using SparseDirectSolver, MadNLP, MadNLPGPU, ExaModels, ExaModelsPower
-using CUDA, CUDSS, Metis, Printf
+using Metis, Printf
+if BACKEND == "cuda"
+    using CUDA, CUDSS
+else
+    using AMDGPU
+end
 
 # the MadNLP linear-solver wrappers (SDSSolver = stock, SDSProtoSolver = PR #107
 # prototype kernels), shipped with the package under bench/e2e
@@ -31,16 +36,22 @@ function solve(label, model, linear_solver)
               fixed_variable_treatment = MadNLP.RelaxBound,
               linear_solver, tol = 1e-6, print_level = MadNLP.ERROR, max_iter = 200)
     madnlp(model; kwargs..., max_iter = 2)                       # warm-up (JIT)
-    GC.gc(); CUDA.reclaim()
+    GC.gc()
+    BACKEND == "cuda" && CUDA.reclaim()
     sol = madnlp(model; kwargs...)
     c = sol.counters
     @printf "%-10s %-20s iters %3d  obj %.8e  wall %6.2f s  linsolve %6.2f s\n" label string(sol.status) c.k sol.objective c.total_time c.linear_solver_time
     return sol
 end
 
-println("device: ", CUDA.name(CUDA.device()))
-model, _ = ac_opf_model("pglib_opf_case78484_epigrids.m"; backend = CUDABackend())
+println("device: ", BACKEND == "cuda" ? CUDA.name(CUDA.device()) : string(AMDGPU.device()))
+model, _ = ac_opf_model("pglib_opf_case78484_epigrids.m";
+                        backend = BACKEND == "cuda" ? CUDABackend() : ROCBackend())
 
-solve("cuDSS", model, MadNLPGPU.CUDSSSolver)
-solve("SDS", model, SDSProtoSolver)
+solvers = BACKEND == "cuda" ?
+    ["cuDSS" => MadNLPGPU.CUDSSSolver, "SDS" => SDSProtoSolver] :
+    ["SDS stock" => SDSSolver, "SDS proto" => SDSProtoSolver]
+for (label, ls) in solvers
+    solve(label, model, ls)
+end
 println("done")

From 23223d8cd09a833ab0be670355e3811a4e0ac438 Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 22:08:06 -0400
Subject: [PATCH 08/16] bench/e2e: fresh_machine.jl develops SDS from an
 editable checkout (~/.julia/dev) instead of a frozen add

---
 bench/e2e/fresh_machine.jl | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/bench/e2e/fresh_machine.jl b/bench/e2e/fresh_machine.jl
index 5ddd50b..e3a603d 100644
--- a/bench/e2e/fresh_machine.jl
+++ b/bench/e2e/fresh_machine.jl
@@ -14,7 +14,11 @@ const BACKEND = get(ENV, "SDS_BACKEND", Sys.which("nvidia-smi") === nothing ? "a
 Pkg.activate(joinpath(homedir(), ".sds-fresh-" * BACKEND))
 Pkg.add(["MadNLP", "MadNLPGPU", "ExaModels", "ExaModelsPower", "Metis", "Printf"])
 Pkg.add(BACKEND == "cuda" ? ["CUDA", "CUDSS"] : ["AMDGPU"])
-Pkg.add(url = "https://github.com/exanauts/SparseDirectSolver.jl", rev = "divfree-chol")
+# SDS as a proper dev checkout (editable, canonical path), not a frozen Pkg.add
+const SDS_DEV = joinpath(homedir(), ".julia", "dev", "SparseDirectSolver")
+isdir(SDS_DEV) ||
+    run(`git clone --branch divfree-chol https://github.com/exanauts/SparseDirectSolver.jl $SDS_DEV`)
+Pkg.develop(path = SDS_DEV)
 Pkg.instantiate()
 
 using SparseDirectSolver, MadNLP, MadNLPGPU, ExaModels, ExaModelsPower

From 3ef916fb53ba5dfb36765f0f4363ee650d9180f9 Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 22:15:09 -0400
Subject: [PATCH 09/16] bench/e2e: fresh_machine.jl runs only the
 SDSProtoSolver workflow

---
 bench/e2e/fresh_machine.jl | 22 ++++++++--------------
 1 file changed, 8 insertions(+), 14 deletions(-)

diff --git a/bench/e2e/fresh_machine.jl b/bench/e2e/fresh_machine.jl
index e3a603d..9728e5d 100644
--- a/bench/e2e/fresh_machine.jl
+++ b/bench/e2e/fresh_machine.jl
@@ -2,18 +2,17 @@
 #
 #   julia +1.13 fresh_machine.jl          (or:  julia +1.13 -i fresh_machine.jl)
 #
-# Installs everything into its own environment (~/.sds-fresh), downloads the
-# pglib case on first use, then solves the 78k-bus ACOPF with MadNLP on the GPU.
-# NVIDIA: cuDSS vs SparseDirectSolver on the same device. AMD: SDS stock vs the
-# prototype kernels (there is no cuDSS; a system ROCm installation is required).
-# The vendor is autodetected; force it with SDS_BACKEND=cuda|amdgpu. The workflow
-# passes when both solvers converge in the same iterations to the same objective.
+# Installs everything into its own environment (~/.sds-fresh-<backend>), downloads
+# the pglib case on first use, then solves the 78k-bus ACOPF with MadNLP on the
+# GPU using SparseDirectSolver (the SDSProtoSolver wrapper) as the linear solver.
+# The vendor is autodetected; force it with SDS_BACKEND=cuda|amdgpu. (AMD needs a
+# system ROCm installation.)
 
 import Pkg
 const BACKEND = get(ENV, "SDS_BACKEND", Sys.which("nvidia-smi") === nothing ? "amdgpu" : "cuda")
 Pkg.activate(joinpath(homedir(), ".sds-fresh-" * BACKEND))
 Pkg.add(["MadNLP", "MadNLPGPU", "ExaModels", "ExaModelsPower", "Metis", "Printf"])
-Pkg.add(BACKEND == "cuda" ? ["CUDA", "CUDSS"] : ["AMDGPU"])
+Pkg.add(BACKEND == "cuda" ? ["CUDA"] : ["AMDGPU"])
 # SDS as a proper dev checkout (editable, canonical path), not a frozen Pkg.add
 const SDS_DEV = joinpath(homedir(), ".julia", "dev", "SparseDirectSolver")
 isdir(SDS_DEV) ||
@@ -24,7 +23,7 @@ Pkg.instantiate()
 using SparseDirectSolver, MadNLP, MadNLPGPU, ExaModels, ExaModelsPower
 using Metis, Printf
 if BACKEND == "cuda"
-    using CUDA, CUDSS
+    using CUDA
 else
     using AMDGPU
 end
@@ -52,10 +51,5 @@ println("device: ", BACKEND == "cuda" ? CUDA.name(CUDA.device()) : string(AMDGPU
 model, _ = ac_opf_model("pglib_opf_case78484_epigrids.m";
                         backend = BACKEND == "cuda" ? CUDABackend() : ROCBackend())
 
-solvers = BACKEND == "cuda" ?
-    ["cuDSS" => MadNLPGPU.CUDSSSolver, "SDS" => SDSProtoSolver] :
-    ["SDS stock" => SDSSolver, "SDS proto" => SDSProtoSolver]
-for (label, ls) in solvers
-    solve(label, model, ls)
-end
+solve("SDS", model, SDSProtoSolver)
 println("done")

From 347c2da5a8f88ca2e407a521b30770bba1de13d5 Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 22:15:43 -0400
Subject: [PATCH 10/16] bench/e2e: fresh_machine.jl ends in a plain madnlp(m)
 call

---
 bench/e2e/fresh_machine.jl | 30 ++++++++++--------------------
 1 file changed, 10 insertions(+), 20 deletions(-)

diff --git a/bench/e2e/fresh_machine.jl b/bench/e2e/fresh_machine.jl
index 9728e5d..de5b15c 100644
--- a/bench/e2e/fresh_machine.jl
+++ b/bench/e2e/fresh_machine.jl
@@ -11,7 +11,7 @@
 import Pkg
 const BACKEND = get(ENV, "SDS_BACKEND", Sys.which("nvidia-smi") === nothing ? "amdgpu" : "cuda")
 Pkg.activate(joinpath(homedir(), ".sds-fresh-" * BACKEND))
-Pkg.add(["MadNLP", "MadNLPGPU", "ExaModels", "ExaModelsPower", "Metis", "Printf"])
+Pkg.add(["MadNLP", "MadNLPGPU", "ExaModels", "ExaModelsPower", "Metis"])
 Pkg.add(BACKEND == "cuda" ? ["CUDA"] : ["AMDGPU"])
 # SDS as a proper dev checkout (editable, canonical path), not a frozen Pkg.add
 const SDS_DEV = joinpath(homedir(), ".julia", "dev", "SparseDirectSolver")
@@ -21,7 +21,7 @@ Pkg.develop(path = SDS_DEV)
 Pkg.instantiate()
 
 using SparseDirectSolver, MadNLP, MadNLPGPU, ExaModels, ExaModelsPower
-using Metis, Printf
+using Metis
 if BACKEND == "cuda"
     using CUDA
 else
@@ -33,23 +33,13 @@ end
 include(joinpath(pkgdir(SparseDirectSolver), "bench", "e2e", "MadNLPSDS.jl"))
 using .MadNLPSDS
 
-function solve(label, model, linear_solver)
-    kwargs = (; kkt_system = MadNLP.SparseCondensedKKTSystem,
-              equality_treatment = MadNLP.RelaxEquality,
-              fixed_variable_treatment = MadNLP.RelaxBound,
-              linear_solver, tol = 1e-6, print_level = MadNLP.ERROR, max_iter = 200)
-    madnlp(model; kwargs..., max_iter = 2)                       # warm-up (JIT)
-    GC.gc()
-    BACKEND == "cuda" && CUDA.reclaim()
-    sol = madnlp(model; kwargs...)
-    c = sol.counters
-    @printf "%-10s %-20s iters %3d  obj %.8e  wall %6.2f s  linsolve %6.2f s\n" label string(sol.status) c.k sol.objective c.total_time c.linear_solver_time
-    return sol
-end
-
 println("device: ", BACKEND == "cuda" ? CUDA.name(CUDA.device()) : string(AMDGPU.device()))
-model, _ = ac_opf_model("pglib_opf_case78484_epigrids.m";
-                        backend = BACKEND == "cuda" ? CUDABackend() : ROCBackend())
 
-solve("SDS", model, SDSProtoSolver)
-println("done")
+m, _ = ac_opf_model("pglib_opf_case78484_epigrids.m";
+                    backend = BACKEND == "cuda" ? CUDABackend() : ROCBackend())
+
+sol = madnlp(m; linear_solver = SDSProtoSolver,
+             kkt_system = MadNLP.SparseCondensedKKTSystem,
+             equality_treatment = MadNLP.RelaxEquality,
+             fixed_variable_treatment = MadNLP.RelaxBound,
+             tol = 1e-6)

From b0b19f2a1852755859106e56fc67c610e2d08647 Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 22:21:47 -0400
Subject: [PATCH 11/16] bench/e2e: drop the redundant madnlp options (MadNLP
 defaults them under SparseCondensedKKTSystem)

---
 bench/e2e/fresh_machine.jl | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/bench/e2e/fresh_machine.jl b/bench/e2e/fresh_machine.jl
index de5b15c..7863b6c 100644
--- a/bench/e2e/fresh_machine.jl
+++ b/bench/e2e/fresh_machine.jl
@@ -38,8 +38,7 @@ println("device: ", BACKEND == "cuda" ? CUDA.name(CUDA.device()) : string(AMDGPU
 m, _ = ac_opf_model("pglib_opf_case78484_epigrids.m";
                     backend = BACKEND == "cuda" ? CUDABackend() : ROCBackend())
 
+# SparseCondensedKKTSystem is required (the SPD Cholesky path); MadNLP defaults
+# equality/fixed-variable treatment correctly for it. tol: condensed default is 1e-4.
 sol = madnlp(m; linear_solver = SDSProtoSolver,
-             kkt_system = MadNLP.SparseCondensedKKTSystem,
-             equality_treatment = MadNLP.RelaxEquality,
-             fixed_variable_treatment = MadNLP.RelaxBound,
-             tol = 1e-6)
+             kkt_system = MadNLP.SparseCondensedKKTSystem, tol = 1e-6)

From 1cdbe094a35c4e5c87f5fa2d264e1e9e12392e52 Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 22:34:58 -0400
Subject: [PATCH 12/16] fix: reapply the kernel optimizations on current main
 without reverting T23

The previous commits copied subtree.jl/front.jl wholesale from a pre-T23
checkout, silently reverting the per-backend local-memory machinery (the
5-entry SUBTREE_LOCAL_SIZES and the matching _with_local_bytes), which made
every analysis fail with 'update _with_local_bytes' on a clean clone. Both
files are rebuilt from origin/main with only the two optimization hunks
applied (division-free scale-first column step; packed-triangle enumeration
of the trailing update).
---
 src/numeric/front.jl   |  2 +-
 src/numeric/subtree.jl | 24 +++++++++++++-----------
 2 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/src/numeric/front.jl b/src/numeric/front.jl
index 1e0008c..672a759 100644
--- a/src/numeric/front.jl
+++ b/src/numeric/front.jl
@@ -55,7 +55,7 @@ end
 end
 
 # column j, second half: rank-1 update of the trailing block with the scaled column
-# (multiply-add only; the divisions live in the first half as one reciprocal)
+# (multiply-add only; the packed-triangle enumeration replaces the guarded square)
 @inline function _front_chol_scale!(L11, st, piv, j, s, li, front_ncols, ::Val{W}, ::Val{WG}) where {W, WG}
     @inbounds begin
         w = front_ncols[s]
diff --git a/src/numeric/subtree.jl b/src/numeric/subtree.jl
index df41008..022425d 100644
--- a/src/numeric/subtree.jl
+++ b/src/numeric/subtree.jl
@@ -139,8 +139,9 @@ end
 end
 
 # column j, first half: check the pivot and scale the column below the diagonal by it
-# (every work item recomputes `d = sqrt(ajj)` from local memory: cheaper than sharing
-# it through a barrier, and it keeps the second half free of divisions)
+# (every work item recomputes `rd` from local memory: cheaper than sharing it through
+# a barrier, and it keeps the second half free of divisions; the diagonal itself is
+# written in the second half, so no work item reads a value another is overwriting)
 @inline function _subtree_chol_update!(buf, ctl, piv, j, li, ::Val{WG}) where {WG}
     @inbounds begin
         if ctl[_ST_STATUS] == 0
@@ -148,10 +149,10 @@ end
             lf = ctl[_ST_LF] - 1
             ajj = real(buf[lf + _packed(j, j, f)])
             if ajj > 0
-                rd = inv(sqrt(ajj))                        # strictly below the diagonal only:
-                for i in (j + li):WG:f                     # the diagonal is replaced in the second
-                    buf[lf + _packed(i, j, f)] *= rd       # half, so no work item reads a value
-                end                                        # another one is overwriting
+                rd = inv(sqrt(ajj))
+                for i in (j + li):WG:f
+                    buf[lf + _packed(i, j, f)] *= rd
+                end
             elseif li == 1
                 ctl[_ST_STATUS] = j % eltype(ctl)          # not positive (or NaN): stop, as LAPACK potrf
             end
@@ -161,8 +162,9 @@ end
 end
 
 # column j, second half: rank-1 update of the trailing front with the scaled column
-# (multiply-add only: the former divide-in-the-update doubled the flop cost in
-# f64 divisions, the bulk of the fused regime-A kernel's time)
+# (multiply-add only: the former divide-in-the-update doubled the flop cost in f64
+# divisions; the packed-triangle enumeration halves the work items of the old
+# guarded full-square loop)
 @inline function _subtree_chol_scale!(buf, ctl, piv, j, li, ::Val{WG}) where {WG}
     @inbounds begin
         if ctl[_ST_STATUS] == 0
@@ -304,10 +306,10 @@ function _launch_subtrees!(N::Numeric{T}, S::Symbolic, nzval, first, count, ::Va
 end
 
 # the local-memory size class as a compile-time constant (explicit branches: no dynamic dispatch, no allocation;
-# every branch is compiled with the first launch, hence the short list of sizes; branch 5 is the error)
+# every branch is compiled with the first launch, hence the short list of sizes; branch 6 is the error)
 @inline function _with_local_bytes(fn, nbytes::Int)
-    length(SUBTREE_LOCAL_SIZES) == 4 || error("update _with_local_bytes")
-    Base.Cartesian.@nif 5 d -> (nbytes == SUBTREE_LOCAL_SIZES[d]) d -> fn(Val(SUBTREE_LOCAL_SIZES[d])) d -> throw(
+    length(SUBTREE_LOCAL_SIZES) == 5 || error("update _with_local_bytes")
+    Base.Cartesian.@nif 6 d -> (nbytes == SUBTREE_LOCAL_SIZES[d]) d -> fn(Val(SUBTREE_LOCAL_SIZES[d])) d -> throw(
         InvalidValueError("no regime-A kernel with $nbytes bytes of local memory; sizes: $SUBTREE_LOCAL_SIZES"))
 end
 

From 4d7e19e648f8718d83eca4d917ef8c1f5a7536aa Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 22:39:06 -0400
Subject: [PATCH 13/16] bench/e2e: fresh_machine.jl loads CUDSS again
 (MadNLPGPU's CUDA extension requires it to activate the GPU condensed-KKT
 path)

---
 bench/e2e/fresh_machine.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/bench/e2e/fresh_machine.jl b/bench/e2e/fresh_machine.jl
index 7863b6c..f3f40b6 100644
--- a/bench/e2e/fresh_machine.jl
+++ b/bench/e2e/fresh_machine.jl
@@ -12,7 +12,7 @@ import Pkg
 const BACKEND = get(ENV, "SDS_BACKEND", Sys.which("nvidia-smi") === nothing ? "amdgpu" : "cuda")
 Pkg.activate(joinpath(homedir(), ".sds-fresh-" * BACKEND))
 Pkg.add(["MadNLP", "MadNLPGPU", "ExaModels", "ExaModelsPower", "Metis"])
-Pkg.add(BACKEND == "cuda" ? ["CUDA"] : ["AMDGPU"])
+Pkg.add(BACKEND == "cuda" ? ["CUDA", "CUDSS"] : ["AMDGPU"])
 # SDS as a proper dev checkout (editable, canonical path), not a frozen Pkg.add
 const SDS_DEV = joinpath(homedir(), ".julia", "dev", "SparseDirectSolver")
 isdir(SDS_DEV) ||
@@ -23,8 +23,8 @@ Pkg.instantiate()
 using SparseDirectSolver, MadNLP, MadNLPGPU, ExaModels, ExaModelsPower
 using Metis
 if BACKEND == "cuda"
-    using CUDA
-else
+    using CUDA, CUDSS      # CUDSS must be LOADED: MadNLPGPU's CUDA extension (the GPU
+else                       # condensed-KKT constructors) activates only with it present
     using AMDGPU
 end
 

From a5be9315266be66edac04a8601a26970bfe9dc87 Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Fri, 9 Oct 2026 22:46:31 -0400
Subject: [PATCH 14/16] bench/e2e: fresh_machine.jl updates an existing dev
 checkout to the branch head

---
 bench/e2e/fresh_machine.jl | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/bench/e2e/fresh_machine.jl b/bench/e2e/fresh_machine.jl
index f3f40b6..0c6f16a 100644
--- a/bench/e2e/fresh_machine.jl
+++ b/bench/e2e/fresh_machine.jl
@@ -13,10 +13,17 @@ const BACKEND = get(ENV, "SDS_BACKEND", Sys.which("nvidia-smi") === nothing ? "a
 Pkg.activate(joinpath(homedir(), ".sds-fresh-" * BACKEND))
 Pkg.add(["MadNLP", "MadNLPGPU", "ExaModels", "ExaModelsPower", "Metis"])
 Pkg.add(BACKEND == "cuda" ? ["CUDA", "CUDSS"] : ["AMDGPU"])
-# SDS as a proper dev checkout (editable, canonical path), not a frozen Pkg.add
+# SDS as a proper dev checkout (editable, canonical path), not a frozen Pkg.add;
+# an existing checkout is updated to the branch head (local edits block the pull loudly)
 const SDS_DEV = joinpath(homedir(), ".julia", "dev", "SparseDirectSolver")
-isdir(SDS_DEV) ||
-    run(`git clone --branch divfree-chol https://github.com/exanauts/SparseDirectSolver.jl $SDS_DEV`)
+const SDS_REF = "divfree-chol"
+if isdir(SDS_DEV)
+    run(`git -C $SDS_DEV fetch origin $SDS_REF`)
+    run(`git -C $SDS_DEV checkout $SDS_REF`)
+    run(`git -C $SDS_DEV merge --ff-only FETCH_HEAD`)
+else
+    run(`git clone --branch $SDS_REF https://github.com/exanauts/SparseDirectSolver.jl $SDS_DEV`)
+end
 Pkg.develop(path = SDS_DEV)
 Pkg.instantiate()
 

From 836ff2496a1f99feaec5fe0fcce7db5d764e6283 Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Sat, 10 Oct 2026 00:04:28 -0400
Subject: [PATCH 15/16] bench/e2e: vendor detection runs the GPU tool
 (nvidia-smi on PATH without a driver chose cuda on AMD nodes); document
 --startup-file=no everywhere

---
 bench/e2e/fresh_machine.jl | 16 ++++++++++++++--
 bench/e2e/fresh_machine.sh |  2 +-
 2 files changed, 15 insertions(+), 3 deletions(-)

diff --git a/bench/e2e/fresh_machine.jl b/bench/e2e/fresh_machine.jl
index 0c6f16a..c9da038 100644
--- a/bench/e2e/fresh_machine.jl
+++ b/bench/e2e/fresh_machine.jl
@@ -1,6 +1,7 @@
 # MadNLP + SparseDirectSolver workflow on a fresh GPU machine, one file:
 #
-#   julia +1.13 fresh_machine.jl          (or:  julia +1.13 -i fresh_machine.jl)
+#   julia +1.13 --startup-file=no fresh_machine.jl
+#   julia +1.13 --startup-file=no -i fresh_machine.jl          (drop into a REPL after)
 #
 # Installs everything into its own environment (~/.sds-fresh-<backend>), downloads
 # the pglib case on first use, then solves the 78k-bus ACOPF with MadNLP on the
@@ -9,7 +10,18 @@
 # system ROCm installation.)
 
 import Pkg
-const BACKEND = get(ENV, "SDS_BACKEND", Sys.which("nvidia-smi") === nothing ? "amdgpu" : "cuda")
+# vendor detection must RUN the tool, not just find it: cluster nodes often have
+# nvidia-smi on PATH without an NVIDIA driver (and vice versa)
+_works(cmd) = try
+    success(pipeline(cmd; stdout = devnull, stderr = devnull))
+catch
+    false
+end
+const BACKEND = get(ENV, "SDS_BACKEND") do
+    _works(`nvidia-smi -L`) ? "cuda" :
+    (_works(`rocm-smi`) || _works(`rocminfo`)) ? "amdgpu" :
+    error("no working GPU tool found (nvidia-smi / rocm-smi); set SDS_BACKEND=cuda|amdgpu")
+end
 Pkg.activate(joinpath(homedir(), ".sds-fresh-" * BACKEND))
 Pkg.add(["MadNLP", "MadNLPGPU", "ExaModels", "ExaModelsPower", "Metis"])
 Pkg.add(BACKEND == "cuda" ? ["CUDA", "CUDSS"] : ["AMDGPU"])
diff --git a/bench/e2e/fresh_machine.sh b/bench/e2e/fresh_machine.sh
index d9d1490..94e329a 100755
--- a/bench/e2e/fresh_machine.sh
+++ b/bench/e2e/fresh_machine.sh
@@ -51,7 +51,7 @@ if ! command -v juliaup >/dev/null 2>&1 && [ ! -x "$HOME/.juliaup/bin/juliaup" ]
 fi
 export PATH="$HOME/.juliaup/bin:$PATH"
 juliaup add 1.13 2>/dev/null || true
-JL="julia +1.13"
+JL="julia +1.13 --startup-file=no"
 $JL --version
 
 # --- 2. Repository ------------------------------------------------------------

From bd9f01f43ae389e2222d44e86d95574fd4ccac71 Mon Sep 17 00:00:00 2001
From: Sungho Shin <sushin@mit.edu>
Date: Sat, 10 Oct 2026 00:06:14 -0400
Subject: [PATCH 16/16] bench/e2e: NVIDIA detection also requires a listed
 device, not just a clean nvidia-smi exit

---
 bench/e2e/fresh_machine.jl | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/bench/e2e/fresh_machine.jl b/bench/e2e/fresh_machine.jl
index c9da038..ef1af65 100644
--- a/bench/e2e/fresh_machine.jl
+++ b/bench/e2e/fresh_machine.jl
@@ -17,8 +17,13 @@ _works(cmd) = try
 catch
     false
 end
+_has_nvidia_gpu() = try
+    occursin("GPU", read(`nvidia-smi -L`, String))   # exit 0 AND at least one device listed
+catch
+    false
+end
 const BACKEND = get(ENV, "SDS_BACKEND") do
-    _works(`nvidia-smi -L`) ? "cuda" :
+    _has_nvidia_gpu() ? "cuda" :
     (_works(`rocm-smi`) || _works(`rocminfo`)) ? "amdgpu" :
     error("no working GPU tool found (nvidia-smi / rocm-smi); set SDS_BACKEND=cuda|amdgpu")
 end