diff --git a/bench/e2e/fresh_machine.jl b/bench/e2e/fresh_machine.jl new file mode 100644 index 0000000..ef1af65 --- /dev/null +++ b/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-), 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) diff --git a/bench/e2e/fresh_machine.sh b/bench/e2e/fresh_machine.sh new file mode 100755 index 0000000..94e329a --- /dev/null +++ b/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 --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" diff --git a/src/numeric/front.jl b/src/numeric/front.jl index 9652732..672a759 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,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 diff --git a/src/numeric/subtree.jl b/src/numeric/subtree.jl index 27031dc..022425d 100644 --- a/src/numeric/subtree.jl +++ b/src/numeric/subtree.jl @@ -138,7 +138,10 @@ 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 @@ -146,17 +149,9 @@ 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 + 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 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