linstabpy v1.0.0rc8

Added --- cuDSS GPU sparse-direct backend

The largest performance change since the package was created. A new GPU
backend factors the linearized operator with NVIDIA cuDSS (through
nvmath-python) and wraps it in the same right-preconditioned FGMRES that
already owns the accuracy of the iterative path. On the 450,000-unknown
SWBLI resolvent operator (2.0--2.3e7 nonzeros, complex double) one
factorization drops from 41--55 s (STRUMPACK on an NVIDIA RTX A6000) to
5.7--9.0 s (cuDSS, double precision) on the same card. A complete
frequency-line resolvent point --- two refactorizations plus ten forward and
ten adjoint solves converged to 1e-10 --- costs 23 s, against 60--82 s for
the four-rank MUMPS CPU line mode and 130--180 s for the STRUMPACK GPU
path. The GPU changes from the slower arm into the fastest solver on the
node, in full double precision, with no accuracy loss.

• Added linalg/cudss_factor.py: a cuDSS solver wrapper with the four
  operations the resolvent map needs: plan (analysis, once per sparsity
  pattern), factorize, refactorize (new values on the kept pattern,
  through reset_operands(a=...) --- nvmath copies a host operand to the
  device when it is set, not on each call), and per-right-hand-side solve
  (reset_operands(b=...) for the same reason). The cuDSS threading layer
  (libcudss_mtlayer_gomp) is discovered inside the installed NVIDIA wheel
  and enabled with the job's host threads; it halves the analysis time.
• Added factor_backend: cudss to the solver configuration: FGMRES in
  complex double, right-preconditioned by the cuDSS LU, tolerance owned by
  the iteration (ksp_rtol, default 1e-10). The direct MUMPS and STRUMPACK
  paths are unchanged.
• Adjoint solves build a second cuDSS factorization of the conjugate
  transpose. cuDSS 0.8 has no transposed or adjoint solve
  (CUDSS_CONFIG_SOLVE_MODE accepts only its default), so the dual factor
  is the design, not a heal: both factors share one analysis cost profile
  and refactorize together along a frequency line.
• Added factor_precision: single (opt-in, cudss only): the factorization
  runs in complex64 --- 0.7--1.2 s per factor, half the device memory ---
  while FGMRES stays in complex double and restores the solution to 1e-12.
  The option force-enables the two settings measured to be necessary:
  MC64-style matching (MAX_DIAG_PRODUCT) and double-precision row/column
  equilibration before the cast. Without matching the single-precision
  factor fails near zero frequency; with both, it converges at every tested
  point of the gain map.
• The adjoint self-test for this backend draws its probe vectors from the
  forcing-masked subspace instead of using dense random vectors. Random
  probes excite the boundary-penalty rows (row norms up to 5.5e12), where a
  single-precision factor cannot precondition; they would fail a
  configuration that is healthy on every right-hand side the resolvent map
  actually produces.
• Added the optional dependency group cudss (pip install "linstabpy[cudss]"): nvmath-python[cu12] >= 1.0, < 2, which brings the
  cuDSS 0.8 wheels. Requires a CUDA 12 runtime and an NVIDIA GPU supported
  by cuDSS. The core install is unchanged.

Performance --- measured (dual-socket Xeon Silver 4210R, MKL, one RTX A6000 48 GB)

All numbers from the 450,000-unknown SWBLI operator of the validation
study, on an otherwise idle node, with the solve tolerance at 1e-10.

• cuDSS double: analysis 3.3--3.4 s (threading layer on; 6.7--8.1 s
  without), factorization 5.7--9.0 s, 4.1--5.8e8 nonzeros in L+U, one solve
  26--30 ms, two resident factors (K and its adjoint) 15--20 GB device
  memory.
• Refactorization on the kept analysis costs a full numeric factorization
  (9.0--9.2 s); what a frequency line saves is the analysis. Timings are
  reproducible across independent runs to about 1 percent.
• Line-mode point (two refactorizations + 20 solves): 23 s double, 12--14 s
  single at a benign point --- the production sweep below shows single is
  slower overall. Reference points on the same operator and node: MUMPS CPU line
  mode 60--82 s (4 ranks), MUMPS factorization alone 24--33 s (10 threads),
  STRUMPACK GPU point 130--180 s.
• FGMRES iteration counts: 1--2 (double), 5--20 (single with matching and
  equilibration, at points where it converges; the 47-point sweep below found
  21 points where it exceeds the 400-iteration cap).
• Production campaign (two A6000s, one job per card, two concurrent): the
  47-point beta=40/45 subset runs at 57.1 s/pt as independent jobs (CPU record
  118.5 s/pt) and 34.0 s/pt in line mode (CPU 76.6 s/pt); the full 200-point
  map runs at 57.1 s/pt, 103.9 min wall against 83.6 min for five concurrent
  CPU jobs on twenty cores --- wall-clock is bound by device count, per-point
  cost halves. Run records committed as
  examples/savarino/savarino_cudss_*.csv.
• Single precision, measured end-to-end on the same 47 points: converged
  points average 207 s/pt against 58 s/pt double (the weaker factor needs the
  extra iterations), and 21 of 47 points --- the resonant band
  2 <= f <= 7 at beta=45, the response peak at beta=40, and the
  zero-frequency corner --- exceed 400 FGMRES iterations. Evaluated and
  rejected for production; the opt-in stays for the masked-RHS niche and as
  the memory route (halved factor storage).

Accuracy --- measured

• Solutions through the cuDSS-preconditioned FGMRES match the MUMPS double
  reference to 1.2e-12 relative; the discrete adjoint identity holds to
  2.5e-14--5.9e-12 with the dual factor.
• The full production campaign reproduces the committed map-of-record gains
  to every printed digit: all 200 map points and both 47-point arms (jobs and
  line mode) show zero difference at the CSV's ten-digit precision.
• The single-precision option reaches the same 1e-12 through the double
  iteration. Its scope is the production right-hand-side family (forcing-
  and response-masked vectors); on dense unmasked vectors the iteration can
  stall, which is why the option stays opt-in and the self-test uses masked
  probes.
• For a unit-norm right-hand side the attainable true-residual floor is
  about 1e-5, set by the boundary-penalty rows; MUMPS shows the same floor
  on the same operator. Solution error stays at 1e-12.

Fixed --- adjoint documentation

• Corrected the explanation of why the GPU path needs a second
  factorization. PETSc's STRUMPACK interface registers no transposed solve,
  and MatSolveTranspose then silently substitutes the forward solve (the
  factor's symmetry flag is unknown, which the fallback treats as true).
  The measured "transpose defect" therefore comes from the PETSc interface,
  not from STRUMPACK's arithmetic. Earlier comments and the paper draft
  attributed it to STRUMPACK's transpose solve; code comments and
  documentation now state the fallback. The dual-factor construction was
  and remains the correct route on both GPU backends.

Fixed --- reproduction documentation

• A reproduction audit of the shipped tree closed four seams between the new
  backend and its documentation: the "Current Status" backend table in
  docs/COMPUTE_BACKENDS.md still said cuDSS was not wired (contradicting
  the new section in the same file); the install line named only the cudss
  extra when the iterative path also needs the petsc extra, and did not say
  that a standard CPU PETSc/SLEPc build is enough (the GPU work happens
  inside cuDSS); the Savarino example README did not mention the two shipped
  cuDSS configs (savarino_point_cudss.yaml, savarino_point_cudss32.yaml)
  or that the campaign generators drive them via --config; and both
  campaign generators' --config help text named only the STRUMPACK GPU
  route. All four are corrected.

Changed --- Savarino base flow fetched, no longer tracked

• The published Savarino base flow (reference/baseflow/*.dat, 10.9 MB) is no
  longer tracked: its Zenodo record (10.5281/zenodo.18854407) is open access
  but states no redistribution license, so the repository ships no file from
  it. fetch_zenodo.py baseflow downloads it (md5 checked against the Zenodo
  files API) and savarino_to_npy.py re-verifies the md5 before converting.
  The shipped reference/savarino_response_*.npz (repackaged from the same
  record) remain tracked pending the same license decision.

Added --- every published number now has its run record and recipe in the tree

A second reproduction audit, this time claim-by-claim against the manuscript,
brought in everything a reader needs to check or rerun the measurements:

• The STRUMPACK GPU campaign records (savarino_strumpack_farm_fp64.csv,
  savarino_strumpack_lines_fp64.csv) and the dedicated CPU 47-point
  references (savarino_cpu_farm_47pt.csv, savarino_cpu_lines_47pt.csv) —
  the previously unshipped halves of the backend comparison.
• The backend benchmarks (examples/savarino/benchmarks/): the five
  measurement CSVs behind the per-backend factorization times and the
  single-precision verdict, together with the harness scripts that produced
  them and the operator-dump tool they consume.
• savarino_campaign_walltimes.csv: wall-clock of every campaign arm from
  Slurm accounting, first job start to last job end, with job-id ranges. This
  corrected one published number: the cuDSS 200-point map wall-clock is
  103.9 min, not 100.2 (the earlier figure was measured from the submission
  loop's completion, while the first jobs were already running).
• savarino_mpi_scaling.csv: MPI strong scaling at the Mack point measured
  on the same node as every other row (190.5 / 154.6 / 148.2 s at 4/8/16
  ranks, gain identical to ten significant digits). This replaces an earlier
  scaling table whose hardware could not be traced to any machine or record.
• compare_backend_modes.py: prints the Chu-weighted correlation between the
  leading response modes of two runs of the same point, so the "GPU modes
  match the CPU modes" statement is checkable with one command.
• The example README now lists, for every shipped CSV, the exact generator
  command and config that produced it.

Evaluated and rejected --- measured dead ends, recorded so they are not re-walked

• Block low-rank (BLR) compression as a loose FGMRES preconditioner, the
  standard advice for cheaper factors: STRUMPACK GPU BLR at tolerances 1e-2
  and 1e-4 produces factors whose plain solve residuals are 1e16--1e31 on
  this operator (the CPU BLR arm reached 4e39) and stalls FGMRES outright,
  while factoring 9--21x slower than the exact GPU factorization (469--1118
  s against 41--55 s). MUMPS BLR at 1e-3 stalls the iteration; at 1e-5 it
  converges in 36--39 iterations and the extra solves cancel the 20 percent
  factor saving. Shock-dominated operators with 11 decades of row-norm
  spread do not compress.
• cuDSS explicit batching (factoring K and its adjoint, or several
  frequencies, in one call): requires device-resident operands and is moot
  at 1-second factors.
• Hybrid CPU+GPU execution and pivot-type overrides: no measurable fit for
  this operator size.
• A single-precision factorization without matching: fails (stalls) near
  zero frequency even with equilibration --- the reason
  factor_precision: single force-enables matching.

Martini compute scripts wrote NPZ and a log only. Plot scripts loaded those
NPZ files and wrote the PNG figures. A conda file installed the CPU stack
from pyproject.toml.

PETSc-free omit-list coverage rose from 75% (4969/6618) to 83% (5467/6618).

Added

• environment.yml installed the package from pyproject.toml on a conda-forge complex PETSc stack.

Changed

• PETSc-only tests skipped collection when mpi4py or petsc4py was missing.
• Auto PETSc matrix layout stayed local when petsc4py was absent.
• driver.run re-exported RUNTIME_NPZ_CONTRACT_KEYS from the payload module.
• Dissipation tests dropped PETSc stubs after import. Later tests saw a missing stack.
• README installation named the CPU stack. It stated that no GPU conda file ships.
• Martini example_*.py drivers wrote NPZ and a log only.
• example_*_plots.py loaded NPZ files and wrote PNG figures.
• run_all_examples.sh computed each Martini case, then ran the plot module.
• benchmark.sh ran Example 0 as examples.martini.example_0_curvilinear_validation.
• README, tutorial, and runner help used examples.martini.* module names.

Removed

• .pre-commit-config.yaml. The file-size guard stayed in pytest and CI.
• Unused tests/fixtures/curvilinear_derivative_parity.json.
• Campaign leftover examples/martini/physics_audit_campaign_vs_repo.json.
• Example 0 write to artifacts/curvilinear_derivative_parity.json.

Verified

• Tests locked environment.yml against pyproject.toml.
• Tests locked the Martini compute/plot split and the leftover-file ban.
• PETSc-free omit-list coverage rose from 75% (4969/6618) to 83% (5467/6618).
• Full-package coverage without PETSc was 50% (5467/11008).
• .coveragerc left out 4390 solver, linalg, and PETSc-helper statements at 0%.
• 31 omit-list files reached 100%. 22 omitted files stayed at 0%.
• The PETSc-free suite passed 1125 tests. It skipped 113 tests. 29 MATLAB-parity tests stayed deselected.
• ruff check and ruff format --check passed on every Python file.
• Tests covered driver CLI, sponge, shock cleanup, Floquet shift, and wavemaker inputs.
• Tests covered grid-field load, filters, conservative-momentum forcing, and typed solver config.
• Tests covered mesh create limits, SBP dissipation errors, and empty Dirichlet recipes.

Omit-list files (.coveragerc source, solvers and linalg omitted). Hit 5467, miss 1151, stmts 6618.

File	Hit	Miss	Stmts	Cover
resolvent/gain.py	19	100	119	16%
operators/spatial.py	39	101	140	28%
stability/_eigenvalues.py	102	232	334	31%
resolvent/operator.py	18	33	51	35%
resolvent/analysis.py	47	78	125	38%
driver/run.py	447	213	660	68%
utils/provenance.py	78	26	104	75%
operators/assembly.py	519	165	684	76%
boundary_conditions/recipe.py	172	44	216	80%
boundary_conditions/dirichlet.py	150	32	182	82%
utils/cost_metrics.py	42	8	50	84%
boundary_conditions/neumann.py	127	24	151	84%
mesh/sbp_generator.py	147	23	170	86%
operators/metrics.py	178	15	193	92%
mesh/mesh.py	615	47	662	93%
driver/io.py	55	2	57	96%
config/typed.py	296	4	300	99%
baseflow/sponge.py	205	2	207	99%
stability/sensitivity.py	104	1	105	99%
baseflow/preprocess.py	139	1	140	99%
__init__.py	28	0	28	100%
baseflow/__init__.py	5	0	5	100%
baseflow/baseflow.py	21	0	21	100%
baseflow/grid_fields.py	63	0	63	100%
baseflow/shock_cleanup.py	215	0	215	100%
baseflow/viscosity.py	23	0	23	100%
boundary_conditions/__init__.py	11	0	11	100%
config/__init__.py	4	0	4	100%
config/defaults.py	3	0	3	100%
config/schema.py	53	0	53	100%
driver/__init__.py	1	0	1	100%
driver/cli.py	63	0	63	100%
driver/payload.py	12	0	12	100%
driver/preflight.py	165	0	165	100%
mesh/__init__.py	11	0	11	100%
mesh/fdweights.py	25	0	25	100%
mesh/filters.py	189	0	189	100%
mesh/sbp.py	211	0	211	100%
metrics_schema.py	38	0	38	100%
mumps_options.py	3	0	3	100%
operators/__init__.py	4	0	4	100%
operators/coefficients.py	453	0	453	100%
operators/energy_norm.py	12	0	12	100%
resolvent/__init__.py	12	0	12	100%
resolvent/conservative_forcing.py	99	0	99	100%
stability/__init__.py	20	0	20	100%
stability/acceptance.py	93	0	93	100%
utils/__init__.py	2	0	2	100%
utils/_solver_cfg.py	67	0	67	100%
utils/floquet.py	32	0	32	100%
utils/yaml_io.py	30	0	30	100%
Total	5467	1151	6618	83%

Files omitted by .coveragerc. Statement count from the coverage parser. Hit 0 without PETSc.

File	Hit	Miss	Stmts	Cover
solvers/router/stable_sweep.py	0	675	675	0%
solvers/router/resolvent_map.py	0	573	573	0%
linalg/factorize.py	0	452	452	0%
utils/petsc_core.py	0	387	387	0%
solvers/eigenvalue.py	0	320	320	0%
solvers/router/map_walk.py	0	265	265	0%
linalg/backends.py	0	262	262	0%
solvers/router/ram_autoselect.py	0	252	252	0%
solvers/router/anchors.py	0	204	204	0%
solvers/router/recycling.py	0	199	199	0%
solvers/router/certification.py	0	172	172	0%
solvers/router/break_even.py	0	159	159	0%
solvers/spatial_eig.py	0	102	102	0%
solvers/router/band_marking.py	0	96	96	0%
solvers/router/map_record.py	0	78	78	0%
utils/petsc_init.py	0	66	66	0%
solvers/_eps_helpers.py	0	63	63	0%
linalg/shifted.py	0	21	21	0%
solvers/router/__init__.py	0	19	19	0%
linalg/__init__.py	0	13	13	0%
solvers/__init__.py	0	11	11	0%
solvers/petsc_utils.py	0	1	1	0%
Total	0	4390	4390	0%