Repository navigation
Performance improvements and fairer comparison - #9
Merged
Merged
Conversation
The power-flow reuse section still described parameter overrides and branch states as the caller's own contract. They are recorded as detached value snapshots at prepare time and compared on every reuse, so a replaced, added, removed or in-place edited entry is rejected instead of silently reusing stale factors, while an equal-valued replacement keeps the factors and its gradients. Move the harmonic preparation argument of assemble_harmonic_system out of the returns block into the parameter list, where the generated API reference shows it as the argument it is rather than as a fourth return value. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
…trategies A caller that keeps a HarmonicFlowSystem between solves has no way to size what it holds against its own memory budget: one entry can be a full [B, H, N, N] matrix plus its factorization plus the key snapshot the validity check compares against. HarmonicFlowSystem.nbytes() walks the retained entries and reports their tensor storage, counting shared storage once, in the same spirit as AnchoredSystem.nbytes(). The new runnable benchmark measures the three solve strategies that decide the cost of a scenario study: what a preparation saves when the network, the operating point or neither change; how keeping one factorization, factoring every scenario's own matrix, solving each assembled matrix without keeping a factor, correcting one factorization with an exact low-rank update and an iterative solve compare at several grid sizes and batch sizes; and how the current-injection fixed point and Newton compare in iterations, wall time and gradient cost at nominal loading and near the fixed point's convergence limit. It probes the batched dense LU once per size, so a LAPACK build that refuses batched complex factorizations above a few hundred rows reduces the comparison instead of aborting it. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
The low-rank device-shunt path decided whether to use a Woodbury correction by looking at the width of an update core it had already assembled. On a scenario batch that core carries the batch, so a feeder with a load on most of its buses built a [B, H, N, N] tensor only to discard it and take the assembled path anyway. A 900-row feeder at 128 scenarios and six orders asked for 9.2 GiB in complex128 before any chunking applied, which is precisely the allocation the chunked assembly exists to prevent; on a 12 GB accelerator the study ran out of memory. The selection rule now reads the touched-row count, which comes from the block row indices alone, and returns without assembling anything when the rule rejects it. The accepted path, the rule itself and every result are unchanged. The benchmark records a strategy a LAPACK build refuses instead of aborting the sweep: batched complex factorizations that need row interchanges fail on some builds, which says nothing about the strategies that still work on that machine. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Comparing every strategy against one reference solution only works if that reference is right. On some LAPACK builds a batched complex LU returns a wrong solution without raising: at 600 and 900 matrix rows with sixteen or more matrices in one call, the batched dense path disagreed with the sparse, low-rank, one-at-a-time and iterative paths by two to three times the solution magnitude, and it was the batched path that was wrong. Every solution is now checked against the matrix it claims to have solved, so a strategy stands or falls on its own residual. Two independent references are formed and the batched dense path is recorded as unreliable for that size when it fails the check, which is also printed with the timings so a run on such a machine cannot be read as if all its numbers were comparable. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
A caller sizing a retained harmonic preparation against its memory budget needs to know that a numerical factor entry holds the matrix twice: once as the value and once as the snapshot its validity check compares against. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Some CPU LAPACK builds get the row interchange of a multi-threaded batched LU wrong. Above roughly 140 rows, factoring or solving a stack of more than one matrix fails with an oneMKL ZLASWP parameter error followed by a rejected factorization, while the same call in one thread is exact. That removes every path that hands LAPACK a stack of matrices rather than one: the operating-point device shunt, a batched voltage node source, a switch-state sweep that assembles each state, and the implicit-function adjoint of a scenario batch. A harmonic study of a 450-row feeder over sixteen scenarios raised on eight threads and converged on one. Those calls now go through one module. Each is measured against a probe right-hand side whose entries are distinct, because a constant one is invariant under any row permutation and would hide exactly this failure. A call that raises or misses is repeated in a single thread, which is remembered for the rest of the process so the discovery is paid once; a single-threaded call that is still wrong raises instead of returning a wrong voltage. The threshold follows the precision that was factored, so the single-precision factors of a mixed-precision solve are judged by their own round-off. A build without the defect keeps its parallel factorization and pays only the check: 11 percent of a batched factorization of sixteen 300-row systems, 6 percent at 450 rows, nothing measurable at 900. A single matrix, the CUDA path and the sparse backend never enter the guard, so the nonlinear solvers' per-iteration back-substitutions are untouched. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
The throughput comparison gave the reference tools eight worker processes and pgml a single thread, on four physical cores rather than eight. It now gives every engine the same eight physical cores and the same eight workers or threads, so the CPU columns differ by engine and not by how the work was spread. pgml's single batched call stays as a second configuration of the same engine, on the performance page rather than in the README figure. The README keeps one figure. It says plainly that pgml is the slowest tool for a few scenarios and that it does not overtake power-grid-model above about a thousand buses, and its caption states the allocation both sides had. The performance page carries the rest: the grid-size, harmonic, memory and cost figures, a crossover table, and a verdict naming where pgml wins and where it loses. Two things a reader should take from it. Most of the batching gain needs no accelerator, and the accelerator costs 2.6 times the eight cores it is compared against, so the cost crossover arrives at a smaller grid than the throughput crossover. And pgml's harmonic batch is bounded by memory on large grids, where OpenDSS is five to nineteen times faster. The renderer draws each figure from the recorded results, rejects a point that was validated less deeply than the data declares, and reads the agreement tolerance from the record rather than assuming the fundamental one. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
The harmonic preparation, the reuse of a factorization and the choice between the fixed point and Newton change the numbers of the comparison, so the report states what each is worth, with the figures behind them. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
A chunked harmonic scenario run built a reusable HarmonicFlowSystem for every grid, so a small feeder paid the preparation's bookkeeping without getting anything back: the value snapshot of the grid and the key comparison per chunk cost more than the one network assembly and factorization they avoid. Measured on one CPU thread at complex128 over orders 1..13 against the uncached call, the preparation runs at 0.85x to 0.98x up to 201 matrix rows, crosses 1.0 between 201 and 255, and reaches 1.10x to 1.18x at 300 rows and 1.28x to 1.64x from 450 rows up. An accelerator gains at every size measured, 1.1x to 2.5x from 99 rows, because the removed work is a larger share of a GPU call. The run therefore prepares on an accelerator always and on the CPU only from `solver.harmonic.preparation_min_rows` rows up, shipped as 256, the measured crossing. Setting it to 0 prepares at every size; setting it above the largest grid in use keeps the uncached path everywhere. Results are unchanged either way: the preparation only decides whether assembly and factors are rebuilt. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
A device shunt on the `operating_point` basis makes Y(h) one matrix per scenario, and retaining its factorization also retains the matrix itself, because the validity check is an exact comparison against a snapshot. That is a second full copy of a [B, H, N, N] system, up to 804 MiB for 32 scenarios of a 300-row feeder at seven orders, plus a clone and a whole-tensor comparison on every call. Only a replay of the identical batch at identical admittances can use the entry, and a stream of scenarios never does. Retention is therefore off by default. Measured against the uncached call on one CPU thread at complex128, the batched operating-point-shunt cases move from 0.82x-0.91x to 0.99x-1.01x. `cache_batched_factors=True` brings the entry back for a caller that really does solve the same batch again, where it is worth 1.3x to 2.6x. Nothing else changes: the factorization is still built, so results and gradients are identical either way. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Above the fundamental a load is a current source in parallel with a shunt admittance built from its power, and that admittance follows the scenario in both engines: OpenDSS rebuilds and refactorizes its system matrix once per scenario and once per order, holding one sparse matrix at a time because it solves the scenarios one after another. The page said only that the two engines do not solve quite the same device model and left the per-order admittance looking like a cost pgml pays alone. The harmonic section now says what both engines build the shunt from, which pgml configuration reproduces OpenDSS and what the alternative costs and buys, with the agreement of each against a live circuit beside its throughput, so a cheaper model is never read as the same answer arriving faster. It also states where the working set is actually bounded: the solver assembles and factors as many scenarios at a time as fit its own budget, so the peak is flat in the batch, and the exact low-rank path is unavailable because a feeder carries a load on half its rows or more. The claim that a 48 GB card limits a study to sixteen scenarios is replaced by what really stops the measured curve, a budget the benchmark applies to the nominal size of the admittance. The preparation figure of the solver-configuration section is redrawn as well. Its axis carries a speed-up, which reads as a multiple, so its ticks now say 1x and 2x instead of 10^0 and 2x10^0 and the decade subdivisions appear as minor grid lines. Its numbers and its caption are unchanged. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
A CPU LAPACK build that gets the row interchange of a multi-threaded batched LU wrong was repaired by re-running the whole factorization in one thread. That is heavier than the defect requires. Setting the thread count also disables the math library's dynamic thread adjustment, which would otherwise drop it to one thread when it is entered from inside another parallel region; the batched factorization then calls the per-matrix routine from several threads at once, each insisting on its full thread count, and above roughly 145 rows the blocked kernel rejects the pivot array it is handed and reports success anyway. A single-matrix factorization is exactly the case the library threads correctly. The guard now escalates: the batched call as asked, then the same stack factored one matrix at a time at the caller's thread count, then the batched call in one thread as the last resort, since the repair leans on the same library. The stage that works is latched for the process and reported by repair_in_use(). Every correctness property is unchanged: the probe right-hand side with distinct entries, the residual threshold keyed to the precision that was factored, and the computation error raised rather than a wrong voltage returned. The repair returns its factors in LAPACK's column-major layout. Stacking them as they come would leave each later back-substitution a transpose to pay, three to four times the cost of the solve itself on a 450-row stack, which is more than the repair saves. lu_solve stays batched and multi-threaded under the per-matrix repair, because the defect is in the factorization's row interchange and a build that also got the back-substitution wrong is caught by the check, which uses it. Sixteen complex128 systems on eight threads: the guarded call takes 10.0 ms at 200 rows, 69.1 at 450 and 274.2 at 900, against 10.0, 102.8 and 707.5 for the single-threaded repair, at the same backward error. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
The ledger entry for the batched dense CPU guard described a single-threaded retry. It now names the mechanism that decides the repair, the two stages a build can be latched into, the column-major layout the repair has to return, and what the check and the repair cost at 200, 450 and 900 rows. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Co-Authored-By: Claude Opus 5 (1M context) <[email protected]> # Conflicts: # assets/PERFORMANCE.md # assets/readme/README.md # assets/readme/provenance.json # assets/readme/solver_strategy_preparation.svg
The memory section now names the three things the table separates. A worker pool holds nine interpreters, each with its own copy of the library, its own model of the grid and its own working arrays, and that is the three gigabytes it stands at before it solves anything. One process holds the same things once and then adds one matrix, one factorisation and one tensor per quantity the iteration carries, each with a scenario axis. The two therefore grow differently with the batch, so the section states the standing cost and the cost per added scenario per tool, and where the two curves cross. The device-memory paragraph is corrected. The part that scales with the square of the grid is the matrix and its factorisation and does not depend on the batch; the part that grows with the batch is linear in the grid, and at four thousand scenarios it is the larger of the two. The closing section gains the reason the large grids go the other way, measured rather than assumed. Above 512 rows the CPU path already factorises sparsely with a fill-reducing ordering that comes within six per cent of the zero-fill optimum on a radial feeder, and assembly and factorisation together are two and a half per cent of a four-thousand-scenario solve. What costs is the work repeated per scenario and per iteration: the back-substitutions, of which a third to a half is moving data between the tensor library and the sparse solver, and the per- iteration passes over a scenarios-by-rows tensor. On the accelerator the factorisation is dense, so each scenario's back-substitution carries the square of the row count, which a batch does not amortise. Four changes would narrow the gap and are listed by what they are worth; a different ordering is not one of them. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Every sparse factorisation scans the whole square of the row count to find the three nonzeros per row, because the admittance is assembled dense and converted afterwards. Shared over a batch of thousands that is two microseconds per scenario; on its own it is almost all of what one factorisation costs. The closing section says so, because it is what makes a matrix that belongs to one scenario rather than to the whole batch expensive. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
…second Every other figure fixes the batch and compares speed. This one fills a stated memory ceiling with scenarios of one grid and records both the batch that fitted and the throughput reached there, because the tools answer that question differently: a per-scenario engine never holds more than the scenario it is working on, while pgml holds the whole batch and gets its speed from it. The measurement is on one workstation with a much smaller accelerator than the rest of the page uses, and it says so, because the ordering it produces is not the ordering of the fixed-batch figures and the reason is the card rather than the metric. Three things follow and are stated: memory is rarely what stops a tool, power-grid-model is first at capacity on every grid measured here, and the metric rewards a property only one of the tools has. The renderer draws the new figure from the published record with the same colour per tool as the others, an open marker where the search reached its own cap or its time limit rather than the memory ceiling. The phase breakdown behind the large-grid diagnosis is published beside it so its numbers can be checked. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
The scipy SuperLU backend received the assembled Y as a dense torch tensor and converted it, so every factorization scanned N^2 entries to find the O(N) that are nonzero and allocated an N x N temporary to do it. On a distribution feeder that scan is 96 per cent of the cost of one factorization at 1,176 rows and 99 per cent at 4,096. A large fundamental batch hides it behind one conversion per batch; a harmonic study whose device shunt follows the operating point cannot, because there the matrix is per scenario and per order. The positions Y can carry are a property of the topology, not of the values: every contribution the assembler scatters is a component-local block, so the union of each branch's two terminal row sets, each appliance's host node rows and the diagonal is a superset of the nonzero positions at any frequency, operating point, branch state and scenario. `pgml.assembly.ybus_structure` builds that pattern from the grid's incidence, and `lu_factor_system` passes it to the sparse backend, which forms each system's compressed-column data array by gathering nnz values against one shared column structure. The fundamental solve and the harmonic orders supply the pattern wherever the backend resolves to sparse; one pattern serves every order, scenario and chunk of a study. Omitting it keeps the dense scan, so the dense, block and CUDA paths, the small-grid path and the backend selection are untouched. Diagonal equilibration leaves the pattern intact and the ideal-slack path restricts it with the same row selection it applies to Y, so the matrix that reaches SuperLU is the same one as before and results and gradients do not move. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
The current-injection fixed point walks its scenarios-by-rows tensor several times per iteration. Three of those walks produce a number the iteration already has: - Both convergence criteria are decided by comparing every row against its threshold and reducing, then the row maximum is formed anyway to report the achieved level. Where the threshold is one number for every row the two are the same statement, so the criterion is read off that maximum instead: `max_i x_i <= t` and `x_i <= t for all i` agree on every input, non-finite rows included. The voltage-update threshold is always one number; the mismatch threshold is, unless the per-row precision floor binds somewhere. - A scenario that met its criteria is held at the iterate that met them, by a select over the whole batch. Until one has finished that select is the identity, and whether one has is already known: the loop reads its own termination flag back every iteration, so both flags now come from that one synchronisation. - The voltage-update norm was reduced over the batch every iteration although only the last iterate's value is reported. A group whose devices are all constant power carries the ZIP law as its identity, so the terminal magnitude, the per-unit ratio and both polynomials are no longer formed for it: `0*r^2 + 0*r + 1` is exactly 1 and `P*1` is exactly `P`, so the currents are unchanged and the evaluation reads the nameplate power directly. The plan resolves this once, from the coefficients. Every scenario's voltages are unchanged bit for bit, and so are the iteration count and the reported mismatch. Measured on a 4096-row radial feeder at complex128 with 4096 scenarios, the device-current evaluation drops by 28 per cent of its own time. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
With the admittance pattern in hand the sparse backend no longer needs the matrix, but two full O(N^2) copies were still made ahead of it: the ideal slack's free-row block, gathered out of Y, and the diagonal equilibration applied to that block. Together they cost an order of magnitude more than the factorization they prepared, which is what kept the factorization phase at 16 per cent of a small-batch solve on a 4,096-row system after the conversion itself was gone. Neither needs the matrix. The symmetric scaling is van der Sluis and reads the diagonal alone, so it is now built from the diagonal (scales_from_diagonal) and applied to the stored entries; the free-row block is a renumbering of the pattern plus a gather of the entries it keeps. One lu_factor_system call on a 4,096-row system falls from 131 ms to 2.0 ms in Norton mode and from 258 ms to 2.2 ms with an ideal slack, and solves the same system to the last bit. What the scaled matrix is still needed for is the mixed-precision residual, the condition estimate and the adjoint of a gradient-carrying solve. It is built eagerly for those and otherwise left to FactoredSystem.materialised_y(), which rebuilds it on demand; a forward solve never asks, so y_mat is a zero-storage stand-in carrying only the shape, dtype and device it reads. A dense factorization now also checks LAPACK's own report. A negative info is a bad call, not a singular matrix: it leaves unusable factors behind, and torch.linalg.lu_factor does not surface it, so the failure appears either inside the next back-substitution or as a converged solve whose voltages solve nothing. It raises instead. A positive info is an exact zero pivot and keeps torch's own error, so a batched harmonic solve still reports a non-finite scenario rather than failing the batch. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
SuperLU reads and writes a column-major `[m, k]` block of right-hand sides. The backend permuted the caller's scenario-major `[k, m]` torch tensor into a row-major `[m, k]` numpy block, scipy transposed that again into its own column-major copy, and the answer was transposed a third time on the way back — three full copies of the block, measured at 30 to 40 per cent of the backend entry on a 1176-row feeder at complex128. None of them is needed, because the scenario-major block a caller already holds IS the column-major block, read transposed. The right-hand sides now keep that layout: the numpy view shares memory with the torch tensor when it arrives contiguous, scipy's own copy becomes a straight one, and a single-factorization solve hands the answer back as a torch view of SuperLU's output buffer. The solution is unchanged bit for bit; the backend entry drops from 1.8 times the bare back-substitution to 1.02 times it. With the layout settled, a single factorization answering many right-hand sides splits them by columns across the worker threads, which is what the class always described and what the batch-shape generalization stopped doing. Chunks are kept at or above the width where a chunk reproduces the undivided call's blocking bit for bit (measured: identical from 128 columns up, different below it), so the split buys 1.9 times at 4096 right-hand sides without moving a digit of the solution. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Co-Authored-By: Claude Opus 5 (1M context) <[email protected]> # Conflicts: # src/pgml/solver/CONTEXT.md # src/pgml/solver/harmonic.py
… not A stack of matrices is verified by its own backward error, so the library's report is read on the single-matrix path, which has no other check. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>
Co-Authored-By: Claude Opus 5 (1M context) <[email protected]> # Conflicts: # src/pgml/solver/CONTEXT.md # src/pgml/solver/harmonic.py # tests/reference/test_sparse_solver.py
Co-Authored-By: Claude Opus 5.5 <[email protected]>
Where the harmonic orders are factored by the CPU sparse backend, each order is assembled as its entries over the topology's sparsity pattern and factored from them, so a study never forms the dense [B, H, N, N] tensor. The stamps accumulate through a pattern accumulator that raises on any entry outside the pattern. The backend is chosen before assembly by the shared decision function; CUDA, the dense and block backends, the low-rank path, mixed precision and callers that ask for the matrix keep the dense route. Solved voltages are bit-identical to the dense assembly at equal chunking, and a float64 gradcheck covers the values path. The chunking budget charges a structural system per pattern entry, and the prepared-system cache keys the structural entries. Co-Authored-By: Claude Opus 5.5 <[email protected]>
Co-Authored-By: Claude Opus 5.5 <[email protected]>
Co-Authored-By: Claude Opus 5.5 <[email protected]>
Replace every workstation and pre-audit number on the performance page and in the README performance section with the campaign measured on one NVIDIA L40S 48 GB and eight physical cores under the audited protocol: every timed call ends with the voltages in host memory, power-grid-model and the pgml single call run in fresh processes, every pool is warmed per worker, the guard re-solves 256 scenarios, and pandapower solves the three-phase CIGRE LV grid. The harmonic comparison uses one ladder and one time budget for every tool. Sentences the new data contradicts are rewritten in both directions, the engine is named by source digest, and the README figures and provenance are regenerated. Co-Authored-By: Claude Opus 5.5 <[email protected]>
Replace the development-machine working-set figure with the controlled comparison of both engine revisions under one harness on one workstation: the host peak of a harmonic study on the 1,176-row network, the memory kill of the previous revision on the 294-row network, and the study rate of the CPU sparse arms at the largest batch both revisions completed, beside the spread of the dense and GPU arms the change does not touch. Co-Authored-By: Claude Opus 5.5 <[email protected]>
The source checksums now match the published result files of the paper branch, and the recorded paper commit contains the harness that measured them. Co-Authored-By: Claude Opus 5.5 <[email protected]>
Co-Authored-By: Claude Opus 5.5 <[email protected]>
Each panel names its grid, and OpenDSS runs the same batch ladder and time budget as every pgml arm instead of two points. The subsection's table, ratios and memory figures now come from that run on the L40S allocation, with pgml's single-process CPU arm beside the GPU. Co-Authored-By: Claude Opus 5.5 <[email protected]>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
A fairer comparison between pgml and all reference libraries includes
sparse harmonic solve and multipe caching tricks improves CPU performance of pgml