Methodology and reproduction

Run from the repository root after following Build, test and document. CPU comparison dependencies use requirements/benchmark.lock; GPU comparisons use requirements/cuda-benchmark.lock in the separate CUDA environment.

Run the comparison suite after installing its dependencies:

.venv/bin/python -m pip install --cache-dir .cache/pip -r requirements/benchmark.lock
make pyprof PYPROF_SUITE=sparse-solvers

Suites and candidates

Suite

Inputs

Default comparison

cpu-dense-f64-v1

Float64 SPD dense: n=64, 256, 1024, 2048

Native / NumPy-out / eager Torch

cpu-csr-f64-v1

2D Poisson CSR: grid=8, 32, 64, 128

Native / SciPy-out / eager Torch

sparse-solvers

Scaled 2D/3D Poisson and one neutral case

Native and official SciPy CG/PCG / SuperLU

sparse-solvers-large

Unscaled 2D grids 64/128 and 3D grids 20/30

Same solver choices; opt-in

cuda-dense-f64-v1 / cuda-csr-f64-v1

The same dense / CSR input fixtures

Explicit native CUDA / cuBLAS / eager Torch CUDA

Dense inputs are symmetric and strictly diagonally dominant, with diagonal range parameter 100. CSR uses the Dirichlet Poisson operator. Scaled sparse cases use A = D P D, with positive diagonal D spanning 1..20. These parameters describe generation, not exact condition numbers. Fixtures have a seeded known solution and zero initial guess; input hashes include the exact matrix/vector bytes, since a seed alone does not fix every library’s random-number behavior.

Implementation candidates share the selected CG recurrence, residual policy and preconditioner. scipy-out is this project’s CG loop using public SciPy CSR matvec and NumPy vector operations; its product allocation/copy is timed. Torch uses eager operations with preallocated prepared workspace. The explicit native-blas candidate uses the private dense BLAS kernel, with no automatic backend crossover.

The strategy suite instead calls official scipy.sparse.linalg.cg and SuperLU splu. Official SciPy has its own recurrence and stopping policy; its iterations can differ from native CG. Candidate suffixes -cg / -pcg select no preconditioner / Jacobi. Direct and iterative answers can have different accepted errors.

Timing boundaries

  • prepared reuses the matrix, output, workspace and Jacobi setup. Each timed solve restores the initial output.

  • api includes output/workspace allocation, native validation and per-call Jacobi setup. Input generation, canonical CSR conversion and library-ready matrix construction stay outside solve timing.

  • CUDA-only session reuses iteration vectors and execution resources, allocating a fresh output and optional Jacobi inverse on every call for native and eager Torch. It does not snapshot/cache the matrix or preconditioner. Output reset and per-call setup are timed; session construction/destruction stay outside sampling. Torch retains its documented caching allocator. Existing api/prepared specifications keep their meanings.

  • Official SciPy CG allocates iterative workspace in both scopes. SuperLU total factors and solves on every call; reuse factors before timing in either scope. CSR-to-CSC conversion is outside both.

  • CUDA starts with device-resident inputs. Completed wall time includes host dispatch, scalar synchronization, output reset and a completion fence. Uploads and correctness downloads stay outside solve timing. This is neither kernel-only timing nor end-to-end host-to-host latency.

CUDA defaults to native execution_mode="host". Select --cuda-execution-mode device explicitly to benchmark native device scalar control; this adds solver.execution_mode="device" to the artifact. Default host specs retain their existing schema and replay behavior. The flag selects only native CUDA execution: the eager Torch reference keeps its own scalar control and records that the native mode does not apply to it. Both modes use the same input fixtures, numerical options, reset and completion boundaries. Device status packets, residual checks and final completion remain timed. Keep distinct mode artifacts and inspect their operation counts before interpreting a latency difference.

Scopes describe different allocation/setup contracts. Keep their artifacts separate: the strict comparator rejects mismatched scopes, even if operation counts match. Session reuse is not an improvement to the fresh API contract. CPU candidates reject session scope. No hidden pool or automatic selection is introduced.

For a session comparison on resident inputs:

make pyprof PRESET=cuda-release PYTHON="$PWD/.venv-cuda/bin/python" \
    PYPROF_SUITE=smoke-cuda-dense PYPROF_SCOPE=session \
    PYPROF_THREADS=1 PYPROF_ROUNDS=30 PYPROF_REPETITIONS=3 \
    PYPROF_ARGS='--device cuda --cuda-execution-mode device --implementations cuda-native,cuda-cublas,torch-cuda-eager --preconditioner jacobi --warmups 3'

Use PYPROF_SCOPE=api or prepared in separate artifacts to compare the other contracts. For CSR, select cuda-native,torch-cuda-eager.

Setup probes and stored-array counts are supplementary measurements. Stored bytes exclude owning snapshots, workspace, allocator and internal solver storage; they are not peak process or GPU memory.

Accuracy and sampling

Default options are standard CG, fused native vector operations, seed 20261015, rtol=1e-8, atol=1e-12, maxiter=2000, and true-residual refresh/restart every 50 updates. The explicit algorithm is a reference formulation with two matvecs per update; it cannot be combined with Jacobi PCG.

Before and after timing, every candidate must converge and independently satisfy ||b - A x||2 <= max(atol, rtol * ||b||2), with only a 64-machine-epsilon relative summation margin. A scaled math.hypot norm guards against overflow/underflow in this check. Invalid answers are excluded and fail the run; missing optional candidates remain visible as ineligible. Matvec counts include initial and final residual evaluations.

New native and same-policy reference results also record restarts: the number of successful search-direction replacements after a true-residual checkpoint. Initial direction setup is not a restart. Convergence, iteration limits or numerical failure before direction replacement do not increment the count; the explicit algorithm reports zero. Counts reset on every solve, including prepared and Session reuse. Before/after gates require matching iterations, matvecs and restart counts, and profiler JSON records all three.

Official SciPy Krylov solvers and direct LU report restarts=null because this native restart definition is unavailable; summaries show n/a. Schema-1 artifacts without the optional field remain readable and are never rewritten or assigned a fabricated zero. Strict block and paired comparisons reject different known counts or known/unknown mixtures. When both sides are unknown, a legacy comparison remains available with restart_parity_verified=false and an explicit limitation: matching iterations/matvecs alone does not establish full numerical-work parity.

Each run starts a fresh process with native thread limits set before numeric imports. Warmups are untimed. A Williams schedule balances candidate positions and predecessor pairs; requested rounds are rounded up to a complete cycle. Each sample times a fixed number of repeated solves. Median and p5–p95 describe round averages, not independent-run confidence intervals or individual-call tail latency.

Use the same hardware, affinity and compiler for comparisons, and inspect actual thread-pool settings. Shared CI checks correctness and collects artifacts; speed thresholds require a stable measurement host.

Run, summarize and replay

make pyprof PYPROF_SUITE=cpu-csr-f64-v1 PYPROF_THREADS=1 \
    PYPROF_SCOPE=prepared PYPROF_ROUNDS=30 PYPROF_REPETITIONS=5 \
    PYPROF_ARGS='--warmups 3'
PYTHONPATH=build/cpu-release/python .venv/bin/python -m cgbench compare \
    build/cpu-release/benchmarks/cpu-csr-f64-v1-t1-prepared.json
PYTHONPATH=build/cpu-release/python .venv/bin/python -m cgbench replay \
    build/cpu-release/benchmarks/cpu-csr-f64-v1-t1-prepared.json \
    --output build/benchmarks/replayed-csr.json

JSON records raw nanoseconds, actual schedules, options, input hashes, residuals/counts, revision, hardware, package versions, threads and the loaded native binary identity. Keep separate process runs separate. Replay verifies regenerated inputs and the saved schedule while recording the current environment; it does not promise identical elapsed times.

Compare two builds

Build the baseline revision and current checkout separately with matching compiler, optimization and backend settings. Generate fresh artifacts for each build with the commands above, alternating their launch order across several process runs. Compare the new files using the staged package:

PYTHONPATH=build/cpu-release/python .venv/bin/python scripts/compare_performance.py \
    --baseline build/benchmarks/baseline-*.json \
    --current build/benchmarks/current-*.json \
    --output build/benchmarks/comparison.json

Its bootstrap resamples artifact blocks and rounds within each block; few process runs provide limited evidence about process variance.

For a fixed eight-process-block comparison with complementary ABBA/BAAB orders, immutable plans and boundary environment records, see Balanced process campaigns. It launches the same public benchmark commands after builds finish and preserves the original 5% gate and all framework controls.

The paired runner generates fresh fixtures directly from a suite. Each case starts two fresh interpreters and keeps them for that case’s adjacent solves, so previous cases cannot change its allocator history. It needs neither previous JSON nor recorded iteration counts.

When comparing native cores across a benchmark-harness change, stage the same Python harness and wrapper for both extensions. This isolates the core; it does not measure changes to the Python wrapper. The following example preserves each extension’s real path for build metadata:

.venv/bin/python - /absolute/baseline/build/cpu-release/python \
    "$PWD/build/cpu-release/python" <<'PY'
import pathlib
import shutil
import sys

current = pathlib.Path(sys.argv[2])
for label, native in zip(("baseline", "current"), sys.argv[1:], strict=True):
    stage = pathlib.Path("build/comparison") / label
    for package in ("cgbench", "cgpy"):
        shutil.copytree(current / package, stage / package, dirs_exist_ok=True,
                        ignore=shutil.ignore_patterns("*.so", "*.pyd", "__pycache__"))
    binary, = pathlib.Path(native).glob("cgpy/_core*.so")
    link = stage / "cgpy" / binary.name
    link.unlink(missing_ok=True)
    link.symlink_to(binary.resolve())
PY

PYTHONPATH=build/cpu-release/python .venv/bin/python scripts/run_paired_performance.py \
    --baseline-pythonpath "$PWD/build/comparison/baseline" \
    --current-pythonpath "$PWD/build/comparison/current" \
    --suite cpu-csr-f64-v1 --scope prepared --implementations native \
    --output build/benchmarks/paired-csr.json

Both builds must match input hashes, solver options and numerical work, and pass the independent residual gate. The default protocol uses five warmup pairs and 50 alternating pairs, batching toward 10 ms per sample with a cap of 1000 solves. It checks the last solve’s status and counts in each timed batch, with residual and input checks before and after sampling. The 95% interval resamples ten contiguous blocks of five pairs. These intervals are conditional on the two worker processes and do not establish variation across process launches or hosts. Per-case coverage is not adjusted for simultaneous comparisons.

Both comparison scripts use a 5% material-regression budget: an interval entirely at or below 1.05 supports no material regression, one entirely above 1.05 indicates regression, and an interval crossing 1.05 is inconclusive. Preserve both fresh-process and paired results when their conclusions differ.

Diagnostic profiling

make pyprof remains the comparison entry point: it launches the shared cgbench harness with independently verified solutions, multiple candidates and recorded sampling schedules. Use scripts/profile_solve.py to inspect one implementation and one exact suite case in a fresh process. It reuses the shared CLI, fixtures, executors and true-residual gates; its samples are diagnostics, not framework comparison scores.

Build the selected preset before profiling, then keep its package path and environment fixed. For a native CPU solve:

taskset -c 0 env PYTHONPATH="$PWD/build/cpu-release/python" \
    .venv/bin/python scripts/profile_solve.py \
    --suite cpu-dense-f64-v1 --case dense-1024 --implementation native \
    --scope prepared --probe solve --threads 1 --warmups 3 --repetitions 10 \
    --output build/profiles/cpu-native-prepared-dense1024.json

--probe solve times the existing adapter’s completed call, including its initial-guess reset. API scope includes per-call result/workspace/preconditioner setup and lifecycle; prepared scope reuses the resources defined by the adapter. Native --probe workspace and --probe prepare instead time public resource construction and destruction, including completion fences. The prepare probe snapshots the operator and sets up the selected preconditioner. They are separate observations: do not subtract independently sampled medians to attribute a solve’s internal costs.

For CUDA, use the project-local CUDA environment and an existing Nsight Systems installation. --cuda-profiler brackets only the repeated diagnostic calls; the script loads the public NVIDIA runtime from that virtual environment and adds an NVTX range per call. For example:

nsys profile --capture-range=cudaProfilerApi --capture-range-end=stop \
    --trace=cuda,cublas,nvtx --sample=none --cpuctxsw=none \
    --output build/profiles/cublas-jacobi-api-dense1024 \
    taskset -c 0 env PYTHONPATH="$PWD/build/cuda-release/python" \
    .venv-cuda/bin/python scripts/profile_solve.py \
    --suite cuda-dense-f64-v1 --case dense-1024 \
    --implementation cuda-cublas --device cuda --cuda-device 0 \
    --scope api --probe solve --preconditioner jacobi --threads 1 \
    --warmups 3 --repetitions 10 --cuda-profiler \
    --output build/profiles/cublas-jacobi-api-dense1024.json

Change the probe to workspace or prepare for a separate resource capture. Fixture creation, input uploads, independent correctness downloads, warmups, metadata collection, NVTX push/pop and profiler start/stop are outside the recorded wall interval. Profiling adds overhead, and CUDA API durations include host waiting; API and kernel totals are not additive components of unprofiled solve latency.

JSON records raw completed-wall samples, input hashes, solver options, before/after correctness gates, adapter setup and environment metadata. Its metric is diagnostic_phase_completed_wall and framework_score is false. For a frozen staged package, optionally pass its identity.json using --stage-identity; loaded module hashes are checked before and after capture. Keep raw JSON and profiler output under build/ and retain unsuccessful captures when reporting evidence.

When investigating host readback, inspect allocation, transfer and completion calls alongside kernels. cudaMemcpyAsync can block for pageable host memory; its name alone does not establish overlap. Preserve stream completion and storage lifetime when comparing readback strategies, as described in the NVIDIA CUDA 12.8.1 synchronization documentation.