Numeric: division-free Cholesky column step in the front kernels

bench/e2e/fresh_machine.jl

@@ -0,0 +1,68 @@
+# MadNLP + SparseDirectSolver workflow on a fresh GPU machine, one file:
+#
+#   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
+# 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
+# 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
+_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
+    _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
+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;
+# an existing checkout is updated to the branch head (local edits block the pull loudly)
+const SDS_DEV = joinpath(homedir(), ".julia", "dev", "SparseDirectSolver")
+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()
+
+using SparseDirectSolver, MadNLP, MadNLPGPU, ExaModels, ExaModelsPower
+using Metis
+if BACKEND == "cuda"
+    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
+
+# 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
+
+println("device: ", BACKEND == "cuda" ? CUDA.name(CUDA.device()) : string(AMDGPU.device()))
+
+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, tol = 1e-6)

bench/e2e/fresh_machine.sh

@@ -0,0 +1,93 @@
+#!/usr/bin/env bash
+# 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.
+#
+# 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)
+#   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, 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
+exec < /dev/null                     # never read the terminal: safe to detach
+
+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
+fi
+export PATH="$HOME/.juliaup/bin:$PATH"
+juliaup add 1.13 2>/dev/null || true
+JL="julia +1.13 --startup-file=no"
+$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 "== $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).
+# 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 via $RUNNER (downloads the pglib case on first use)"
+$JL --project=bench/e2e -e "using Pkg; Pkg.develop(path = \"$SDS_DIR\"); Pkg.instantiate()"
+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"

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,20 @@ 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 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]
         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) ÷ 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
-            li == 1 && (L11[_packed(j, j, W)] = d)
         end
     end
     return nothing

src/numeric/subtree.jl

@@ -138,25 +138,20 @@ 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 `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
             f = ctl[_ST_F]
             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
+                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
@@ -166,17 +161,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 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
             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) ÷ 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
-            li == 1 && (buf[lf + _packed(j, j, f)] = d)
         end
     end
     return nothing

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