Module bde 

Block-diagonal embedding: solve B scenarios as one sparse system.

Every scenario in a batch shares one topology, so every scenario’s Jacobian shares one sparsity pattern. Rather than B separate sparse solves, stack them:

        ┌ J₁              ┐   ┌ Δx₁ ┐   ┌ f₁ ┐
        │    J₂           │   │ Δx₂ │   │ f₂ │
  J  =  │       ⋱         │ , │  ⋮  │ = │ ⋮  │
        └             J_B ┘   └ Δx_B┘   └ f_B┘

This is the design a GPU path would adopt, and the reason it matters is that it needs no batched solver API. You hand the library one ordinary sparse matrix. That is what makes the AMD path viable even though rocSOLVER’s csrrf_* refactorization routines are not batched, and it is why the same code can later target cuDSS, rocSOLVER, or anything else that factors a sparse matrix.

The equivalence is exact, not approximate:

• No fill crosses a block. Fill at (i, j) needs a path i → k → j through eliminated vertices; the graph is B disconnected components, so no such path exists between blocks. Total fill is the sum of per-block fills.
• Pivoting stays inside its block. A column in block s has nonzeros only in block s’s rows, so partial-pivot search can never reach another scenario.
• The elimination tree is a forest of B independent trees, which is precisely the wide parallelism a GPU refactorization needs — and the reason a single 18k-unknown Jacobian cannot fill a GPU but 256 stacked ones can.

This module exists to prove all of that on real sparse code, on the CPU, before any GPU spend. scripts/bench/jax_oracle.py already established it with dense linear algebra; this extends the result to the actual sparse backends the GPU path will mirror.

It is deliberately not faster than crate::batch::BatchSolver on a CPU — one big factorization beats B small ones only when the hardware wants wide independent work, which is a GPU property. Treat this as an architecture validator and the host-side half of a GPU path.

Structs

BdeScenarioResult

One scenario’s outcome from solve_batch_block_diagonal.

BlockDiagonal

The stacked sparsity pattern for n_scenarios copies of one topology’s Jacobian, plus the per-block recipe needed to refill it.

Functions

solve_batch_block_diagonal

Runs Newton-Raphson on all scenarios at once through a single stacked sparse factorization per iteration, with per-scenario convergence masking.