Sparse and GPU numerics

The linear algebra under the solvers lives in src/GridSim.Core/Numerics/ (CSR storage, sparse LU, ILU(0), BiCGSTAB) and src/GridSim.Gpu/ (the ROCm, OpenCL and CUDA backends behind the ILinearSolver/ISparseLinearSolver seams). This chapter describes each structure and algorithm as implemented, then the measured reason the production power-flow path stays on the CPU.

CSR structures

`CsrMatrix` - real sparse matrices

src/GridSim.Core/Numerics/Sparse.cs:CsrMatrix. Compressed sparse row: RowPtr[n+1], ColIdx[nnz], Val[nnz], columns sorted within each row, plus a cached DiagPtr[n] (index of each row's diagonal entry, -1 if absent). This is what the sparse Newton Jacobian and the estimator's gain matrix are assembled into: each row touches only a bus's electrical neighbours, so the dense O(n^2) storage collapses to O(nnz).

CsrMatrix.Builder accumulates (row, col, value) triplets (zeros dropped at Add), then Build(ensureDiagonal = true):

counting sort into rows, per-row sort by column;
duplicates summed - so stamping is pure accumulation, exactly like the dense +=;
optionally injects an explicit zero diagonal entry into any row lacking one, so DiagPtr is always usable by ILU(0) and the Jacobi preconditioner.

Multiply(x, y) is the straightforward CSR SpMV used by BiCGSTAB.

`ComplexCsr` - the sparse Y-bus

Same pattern with two value arrays G[]/B[] sharing one structure - the network admittance matrix in sparse form. The sparse power flow walks a bus's Y-bus row to compute injections and to emit Jacobian entries, so per-iteration cost is O(nnz), not O(n^2).

The direct sparse LU (Gilbert-Peierls)

src/GridSim.Core/Numerics/SparseLu.cs + src/GridSim.Core/Numerics/SparseLuFactorization.cs. A sparse direct solver: left-looking Gilbert-Peierls LU with partial pivoting over a fill-reducing reverse Cuthill-McKee (RCM) column ordering. Unlike an iterative method it solves Jx = b exactly in one pass regardless of how stiff the Jacobian is (short cables -> huge admittances -> ill-conditioning), which is what converges the full node-level GB replica.

The factorisation is L*U = P*A*Q:

Q is the RCM column permutation, P the partial-pivoting row permutation (stored as pinv: original row -> pivot step);
L is unit lower-triangular (unit diagonal stored as the first entry of each column), U upper-triangular (pivot stored last);
both factors are stored CSC, with capacity grown by doubling as fill appears.

Per column k (the left-looking step): the new column is x = L \ A(:, q[k]) - a single sparse triangular solve whose non-zero pattern is discovered first by a depth-first reach over the graph of the already-built L (SparseLuFactorization.Reach/Dfs, iterative, emitting nodes in topological order), so the numeric work is proportional to the true post-fill non-zeros. The pivot is the largest |x[i]| over not-yet-pivotal rows; a pivot below 1e-300 (or no candidate) returns null - the caller's "singular" signal. Entries above the pivot go to U, the rest are scaled by the pivot into L. Solve(b) then permutes by P, forward-solves L, back-solves U, and un-permutes by Q. (The comments note this replaced an earlier Dictionary<int,double>-row scheme that paid hash overhead and re-derived the sparsity pattern on every solve.)

`Factor` / `Ordering` / `FactorPrepared` - the reuse tiers

SparseLu exposes three entry points, each dropping repeated work for callers with more structure:

Entry point	What it does	Who uses it
`Solve(a, b)`	Factor + one solve.	`CpuSparseLuSolver` - each sparse Newton step (the Jacobian's values and pattern change as Q-limits switch buses).
`Factor(a)`	CSR->CSC conversion (explicit zeros dropped) + RCM + numeric factorisation, returned for repeated right-hand sides.	The estimator's residual-covariance path - one gain matrix solved against one column per measurement.
`Ordering(a)` then `FactorPrepared(n, ap, ai, ax, q)`	RCM alone, then numeric-only refactorisation of an already-prepared CSC matrix with that precomputed ordering.	`SparseWlsStep` - the gain matrix's sparsity pattern is structural (the measurement topology), fixed for a whole bad-data pass, so `Begin` computes the CSC structure, per-row value slots and the ordering once, and each Gauss-Newton iteration only refills values and refactors.

The RCM story

SparseLu.ReverseCuthillMcKee orders the (structurally symmetrised, deduplicated) matrix graph: components are started from minimum-degree nodes, each traversed breadth-first with a node's unvisited neighbours enqueued in ascending-degree order, and the final order reversed. RCM concentrates the band, which for grid-shaped graphs keeps LU fill small; it is computed from the structure only, which is exactly why SparseWlsStep can hoist it out of the iteration loop while the power flow recomputes it per step.

ILU(0) and BiCGSTAB - the iterative fallback

src/GridSim.Core/Numerics/Sparse.cs:Ilu0 is an incomplete LU that keeps A's sparsity pattern (no fill): the standard IKJ sweep over each row's strictly-lower entries, with a tiny-pivot guard (a diagonal below 1e-12 in magnitude is replaced by +/-1e-12). Apply performs the unit-lower forward solve and upper back solve in place.

src/GridSim.Core/Numerics/Sparse.cs:BiCgStab is preconditioned BiCGSTAB for the nonsymmetric systems a Newton step produces: only SpMVs, dots and axpys, so it scales where a dense factorisation would not fit. It reports failure (false) on breakdown (|rho| < 1e-300), stagnation (omega = 0) or iteration exhaustion; convergence is ||r||/||b|| < tol.

In the sparse Newton solver these are strictly the fallback: the direct LU is tried first, and only a null (singular) result routes the step to ILU(0)/BiCGSTAB with tol = min(1e-10, Tolerance*1e-3) and maxIter = 4m + 200 (src/GridSim.Core/PowerFlow/SparseNewtonRaphsonPowerFlow.cs:SparseNewtonRaphsonPowerFlow.RunNewton).

The GPU backends

GridSim.Gpu supplies device implementations of both linear-solver seams across three GPU backends - ROCm/HIP (native AMD Instinct), OpenCL (portable open-driver) and CUDA (NVIDIA) - the first two AMD-first by design. OpenCL and CUDA share one code path: the same ILGPU kernels compiled to either accelerator. Two invariants hold across every GPU solver:

Probing never throws. GpuContext (ILGPU, driving CUDA + OpenCL) and RocmContext (HIP) lazily probe once; on any failure - no device, no driver library, an ILGPU exception, or an OpenCL device without FP64 - they mark the backend unavailable and every solver transparently defers to its CPU counterpart, so the types are always safe to construct and call on any machine.
All device access is serialised through a static Lock. There is one GPU in the box and the CLI fans solves out over Parallel.ForEach, so every allocation, copy, kernel launch and synchronise takes the lock.

ILGPU / CUDA + OpenCL

src/GridSim.Gpu/GpuContext.cs:GpuContext - owns the single ILGPU Context and an Accelerator per accelerator type (CUDA and OpenCL), keyed via AcceleratorFor/IsAvailableFor/DescriptionFor(type) over one shared Lock. Its diagnostics include a direct CUDA-driver P/Invoke probe (cuInit/cuDeviceGetCount, independent of ILGPU) and Fp64Note, a heuristic naming whether the device runs FP64 at data-centre rate or the consumer 1/32-1/64 throttle. OpenCL FP64 is optional (cl_khr_fp64), so an OpenCL device is only accepted after a tiny on-device smoke kernel (1.5 + 1.5 == 3.0) proves FP64 actually compiles and executes - a non-FP64 iGPU is rejected rather than selected-then-faulted.
The two ILGPU solvers below take an AcceleratorType ctor argument (CUDA by default, or OpenCL); the same kernels drive either accelerator, so OpenCL adds a portable AMD/Intel/NVIDIA path with no new kernel code.
src/GridSim.Gpu/GpuDenseLinearSolver.cs:GpuDenseLinearSolver (ILinearSolver) - Gaussian elimination with partial pivoting over an augmented row-major [A | b] on the device. The pivot search is a one-block shared-memory argmax reduction (256 threads) whose result is consumed on-device by the swap kernel - the pivot index never touches the host. The whole forward elimination is queued on one stream (kernels execute in submission order), so the host synchronises exactly twice: after elimination (to read the device-side singularity flag, which preserves the CPU contract of throwing on a singular matrix) and after back-substitution (to read x).
src/GridSim.Gpu/GpuSparseLinearSolver.cs:GpuSparseLinearSolver (ISparseLinearSolver) - Jacobi (diagonal) preconditioned BiCGSTAB from ILGPU kernels: CSR SpMV, the vector updates, and one-block dot-product reductions that return a single scalar to the host (the only points BiCGSTAB genuinely needs host values). Breakdown, stagnation, iteration exhaustion or a non-finite result all return null, so the sparse Newton caller's fall-back-to-direct-CPU contract is identical to the CPU path. A missing/zero diagonal bails to CpuSparseLuSolver up front.

Batched dense solve

IBatchedLinearSolver (src/GridSim.Core/Abstractions) is a second dense seam for the case the single-solve path handles badly: many independent, same-sized systems. bool[] SolveBatch(int count, int n, double[] aFlat, double[] bFlat) solves count n x n systems in one call - aFlat/bFlat are the row-major matrices and right-hand sides concatenated end to end, each system's solution overwriting its own bFlat segment. The return is a per-system success mask: a singular system's segment is filled with NaN and its flag set false, but the batch as a whole never fails, so one bad contingency does not sink the sweep.

src/GridSim.Core/PowerFlow/CpuBatchedLuSolver.cs:CpuBatchedLuSolver - the reference/fallback: loops DenseLuSolver across the batch in parallel. This is the ground truth the GPU path is validated against.
src/GridSim.Gpu/GpuBatchedDenseSolver.cs:GpuBatchedDenseSolver - the same Gaussian-elimination-with-partial-pivoting kernels as GpuDenseLinearSolver, with a batch index threaded through every kernel so one flat augmented buffer carries the whole batch. The decisive difference is amortisation: the entire batch is factored under one GpuContext.Lock acquisition and one pair of host syncs, so the per-solve kernel-launch and global-lock overhead - the cost that makes back-to-back single solves slow - is paid once across count systems rather than once each. With no device present it falls back to CpuBatchedLuSolver, so it is always safe to construct; ROCm has no native batched entry point yet, so ROCm also uses the CPU batched reference.

The motivation is the contingency sweep: a set of same-sized outages collapses from hundreds of individually locked single solves into a handful of batched dispatches. Batching a full AC N-1 across NR iterations that diverge in count, and memory-chunking batches too large for device memory, are future work.

ROCm / HIP

src/GridSim.Gpu/RocmContext.cs:RocmContext - one cheap hipGetDeviceCount P/Invoke inside a broad catch; on a machine without ROCm the first call throws DllNotFoundException and the backend reports unavailable without ever loading native libraries (which is why the assembly builds and runs on a non-AMD box). On a real device it also binds it (hipSetDevice) and reads the true device name via the struct-free hipDeviceGetName (deliberately avoiding the version-sensitive hipDeviceProp_t struct), with an Instinct-aware FP64 note.
src/GridSim.Gpu/RocmDenseLinearSolver.cs:RocmDenseLinearSolver - LU via the native hipSOLVER API (hipsolverDgetrf/Dgetrs, explicit work/lwork, not the ambiguous cuSOLVER-compat hipsolverDn* aliases). A is transpose-packed into column-major (LAPACK layout); the workspace is sized as the max of both buffer-size queries; devInfo is checked and a singular factorisation throws, matching the CPU solver.
src/GridSim.Gpu/RocmSparseLinearSolver.cs:RocmSparseLinearSolver - Jacobi preconditioned BiCGSTAB run fully device-resident: the CSR matrix and every Krylov vector stay on the GPU, the O(nnz) SpMV runs via rocsparse_dcsrmv, the vector algebra via rocBLAS-1 (ddot/daxpy/dscal/dcopy), and the Jacobi apply is an on-device Hadamard product via rocblas_ddgmm. Only the convergence scalars and the final solution cross the bus - no per-iteration vector round-trips.

The native ROCm P/Invoke signatures are matched against the ROCm 6.x headers (each entry carries a // verified vs ROCm <ver> citation); the path is complete and covered by hardware-free seam + ABI-shape tests, awaiting its first runtime run on an AMD Instinct box (via gpucheck + gpu-bench.sh).

Selection

src/GridSim.Cli/SolverBackend.cs:SolverBackend.Select maps --cpu (default), --gpu (auto, AMD-first: ROCm, then OpenCL, then CUDA, then CPU) or --gpu=rocm|opencl|cuda onto the seams, warning and falling back to CPU when the requested accelerator is absent (or an OpenCL device lacks FP64). The AMD-first / open-driver-before-proprietary ordering follows from the FP64 rationale below. SolverBackend.DenseBatch() is the parallel entry point for the batched seam above - it returns a GpuBatchedDenseSolver for a CUDA/OpenCL backend, else the CpuBatchedLuSolver. Nothing outside the CLI ever constructs a GPU solver.

Why power flow stays on the CPU

The CPU backends are not a stopgap while a GPU port matures; for steady-state power flow the CPU is the right processor by the shape of the problem, and the repo carries the measurements (docs/benchmarking.md sec. 4, docs/server-results/). The technical core of the argument:

The matrix is almost entirely zeros, and gets sparser with scale. A bus connects to a physical handful of neighbours regardless of network size, so nnz(Y) ~ n + 2*branches grows linearly while the dense footprint grows quadratically. For gb-full (3,539 buses / 5,918 branches) that is ~15.4 k non-zeros against 12.5 M dense entries - ~ 99.9 % sparse - and a 10,000-bus case is sparser still. The Jacobian inherits the same graph (four blocks interleaved), so its density tracks Y's.

The gold-standard solve is the worst possible GPU kernel. The exact, stiffness-proof method is the sparse-direct LU above - and it is sequential by data dependency (each column's elimination depends on previous fill; the elimination tree for a grid-shaped matrix is deep and narrow), branchy and latency-bound (data-dependent pivoting and fill over irregular index structures), and arithmetically tiny (a whole factorisation is a few million FLOPs). A GPU is built for the opposite: dense, regular, massively parallel FP64 streams. On top of that, Newton is itself iterative: a fresh Jacobian must cross the PCIe bus and a kernel launch must be paid per iteration.

The measurements confirm it at every scale. Summarising docs/server-results/a100-results.md (each GPU compared against its own co-located EPYC host, uncontended, on a validated driver stack):

cpu-sparse beat every GPU engine on every case by one to three orders of magnitude - on case2869pegase, 68 ms (cpu-sparse) vs 4,851 ms (cuda-sparse) vs 54,886 ms (cuda-dense, itself slower than the CPU dense at 46,695 ms).
The H100, with ~3.5x the A100's FP64 throughput, was only ~1.1-1.2x faster than the A100 on the same work - the signature of a launch/transfer-bound workload, not an FP64-bound one. More arithmetic lanes cannot shorten a serial, latency-dominated critical path.
Batched N-1 screening - the one embarrassingly-parallel power-flow workload - also went to the CPU: 28 EPYC cores each running the fast serial solver screened case118 at 932 solves/s against the GPU's 18 (~52x), because independent contingencies fan out across CPU threads (ContingencyAnalysis.ScreenBranchesParallel, safe because GridModel is immutable) while the GPU runs slow solves back to back.
An earlier A6000 run (consumer FP64, 1:64 throttle) had the dense GPU solve ~500x slower than the CPU on a 300-bus case; the A100/H100 runs exist precisely to show that removing the throttle does not change the verdict.

Correctness was never the issue. The GPU backends' converged voltages are bit-for-bit identical to the CPU's on every case where both run, with matching iteration counts. That makes a GPU useful as an independent cross-check - a second implementation arriving at the same answer - but not as a performance path.

The FP64 rationale and the AMD-first ordering. For the parts of the workload that do have dense arithmetic, FP64 rate is the deciding spec. Consumer/workstation NVIDIA parts throttle FP64 (1:32-1:64); AMD Instinct parts run it near full rate. Hence SolverBackend's auto order: ROCm first, then the open-driver OpenCL path, CUDA last, CPU always the default. But the ordering only optimises the case where a GPU is used at all - even an MI300-class card cannot beat one CPU core on a matrix that factors in tens of milliseconds, because the workload has almost no parallel arithmetic to feed it. Every mature power-flow tool (MATPOWER, PSS/E, PowerWorld, pandapower) runs its power flow on the CPU for the same reason. Sizing hardware for a power-flow workload means buying CPU cores and memory, not accelerators.