diff --git a/CHANGELOG.md b/CHANGELOG.md index 795b958e..e3e0911c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,12 +7,14 @@ and this project follows [Semantic Versioning](https://semver.org/). ## [Unreleased] -- Long weighted Gram diagonals and cross-table cells use compensated accumulation, avoiding roundoff-driven Laplacian validation failures and loss of signed cell contributions. +- Long weighted Gram diagonals use compensated accumulation. ### Changed - **BREAKING:** `schwarz_precond::PartitionWeights::Uniform` is a unit variant; `PartitionWeights::len` and `is_empty` are removed. Uniform weights take their shape from `SubdomainCore::n_local`; only explicit weights carry a length. The serialized uniform variant no longer carries a count; the wire-format fixture moves v20 → v21; v20, v19 and v12 remain rejected-payload controls. The fixture number is not a runtime serialization envelope (#425). - +- Schwarz subdomains are factored by approx-chol 0.7 from each channel pair's factor-width-2 Gram factor `V`, one column per observation: sign folding, dominance scaling, floating detection and the Schur complement move into approx-chol, and each irreducible block is a subdomain (#481). +- **BREAKING:** `ApproxSchurConfig::seed` (Python `seed=`) is removed; the Schur sampler draws from `ApproxCholConfig::seed`. +- **BREAKING:** `ScalingConfig::tolerance` must lie in `[0, 1)`, else `BuildError::InvalidScalingTolerance`. - Levels with no cross-factor entries in a channel pair are batched into diagonal Schwarz subdomains per side instead of one subdomain each, cutting per-subdomain build and apply overhead without changing the preconditioner. - **BREAKING:** `ReductionStrategy` is removed: the `reduction` field of `PreconditionerConfig::Additive` and `::Adaptive`, Python's `reduction=` argument and `within.config.ReductionStrategy`, and the strategy argument of `SchwarzPreconditioner::new`/`with_n_dofs` with its `reduction_strategy()` getter. Overlapping subdomain outputs always sum in subdomain order (#472). - **BREAKING:** `schwarz_precond::LocalSolver::solve_local(x, workspace, allow_inner_parallelism)` solves in place, and the required `workspace_len()` replaces `scratch_size()`. `SubdomainEntry::{core, is_empty, scratch_size, apply_weighted_into_with_scratch, apply_weighted_into_atomic}`, `SubdomainCore::{restrict_weighted, prolongate_weighted_add, prolongate_weighted_add_atomic}`, `BuildError::ScratchSizeTooSmall` and `SolveError::Synchronization` are removed. diff --git a/Cargo.lock b/Cargo.lock index 2446e90c..1ca946bb 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -25,12 +25,13 @@ checksum = "940b3a0ca603d1eade50a4846a2afffd5ef57a9feac2c0e2ec2e14f9ead76000" [[package]] name = "approx-chol" -version = "0.6.0" +version = "0.7.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e3f7e4049cc6e5824e9c88f88f98f3759baee45dc5e3fe33fc2e99889414755b" +checksum = "b20c64c7b87a2a7518f8ac03b9797795d9408e33a7d081528985131636d64776" dependencies = [ "num-traits", "rand 0.10.3", + "rayon", "serde", ] diff --git a/README.md b/README.md index 9c6d2037..5a14330f 100644 --- a/README.md +++ b/README.md @@ -161,7 +161,7 @@ as `.build_duration_seconds`; both are preserved across pickling and reuse. | `LocalSolverConfig(approx_chol?, schur?, dense_threshold=24, scaling?)` | Schur reduction + approximate Cholesky. Omit `schur` for the library-default approximate variant; pass `schur=Schur.exact()` to request an exact Schur (slower, used for validation). | | `Schur.approximate(config?)` / `Schur.exact()` | Schur-reduction mode passed as `LocalSolverConfig(schur=...)`. | | `ApproxCholConfig(seed=0, split_merge=None)` | Approximate Cholesky parameters. | -| `ApproxSchurConfig(seed=0, split=1)` | Approximate Schur complement sampling parameters. | +| `ApproxSchurConfig(split=1)` | Approximate Schur complement sampling parameters; draws from `ApproxCholConfig.seed`. | ### Result types @@ -284,7 +284,7 @@ across solves like any other. |---|---|---| | `within::config` | public | `LsmrOptions`, `PreconditionerConfig`, `LocalSolverConfig`, `SchurMode`, `ApproxCholConfig`, `ApproxSchurConfig`, `ScalingConfig` | | `within::error` | public | `WithinError`, `BuildError`, `SolveError` | -| `block_elim` / `channel` / `csr_block` / `domain` / `linalg` / `operator` / `solver` | `pub(crate)` | implementation layers — public items are re-exported at the crate root | +| `block_solver` / `channel` / `domain` / `linalg` / `operator` / `solver` | `pub(crate)` | implementation layers — public items are re-exported at the crate root | ## Project structure diff --git a/benchmarks/_framework.py b/benchmarks/_framework.py index 3e1548a4..70ea2466 100644 --- a/benchmarks/_framework.py +++ b/benchmarks/_framework.py @@ -207,13 +207,13 @@ def standard_solver_configs( """ schur = LocalSolverConfig( approx_chol=ApproxCholConfig(seed=opts.seed), - schur=Schur.approximate(ApproxSchurConfig(seed=opts.seed)), + schur=Schur.approximate(ApproxSchurConfig()), ) # "2-2": both densification knobs on — split_merge=2 (AC2 reduced-system # factorization) and split=2 (denser Schur clique-tree sampling). schur_2_2 = LocalSolverConfig( approx_chol=ApproxCholConfig(seed=opts.seed, split_merge=2), - schur=Schur.approximate(ApproxSchurConfig(seed=opts.seed, split=2)), + schur=Schur.approximate(ApproxSchurConfig(split=2)), ) return [ SolverConfig( diff --git a/benchmarks/suites/ac_comparison.py b/benchmarks/suites/ac_comparison.py index 996f3af2..8eb70dbc 100644 --- a/benchmarks/suites/ac_comparison.py +++ b/benchmarks/suites/ac_comparison.py @@ -36,7 +36,7 @@ def _schur(seed: int, split_merge: int | None) -> LocalSolverConfig: return LocalSolverConfig( approx_chol=ApproxCholConfig(seed=seed, split_merge=split_merge), - schur=Schur.approximate(ApproxSchurConfig(seed=seed)), + schur=Schur.approximate(ApproxSchurConfig()), ) diff --git a/benchmarks/suites/fixest_comparison.py b/benchmarks/suites/fixest_comparison.py index 9cbee70c..fe5a8a85 100644 --- a/benchmarks/suites/fixest_comparison.py +++ b/benchmarks/suites/fixest_comparison.py @@ -57,13 +57,13 @@ def run_fixest_comparison(opts: SuiteOptions) -> list[BenchmarkResult]: ) approx_schur = LocalSolverConfig( approx_chol=ApproxCholConfig(seed=0, split_merge=8), - schur=Schur.approximate(ApproxSchurConfig(seed=0)), + schur=Schur.approximate(ApproxSchurConfig()), ) # "2-2": both densification knobs on — split_merge=2 (AC2 factorization) # and split=2 (denser Schur clique-tree sampling). schur_2_2 = LocalSolverConfig( approx_chol=ApproxCholConfig(seed=0, split_merge=2), - schur=Schur.approximate(ApproxSchurConfig(seed=0, split=2)), + schur=Schur.approximate(ApproxSchurConfig(split=2)), ) solver_configs = [ diff --git a/crates/within-py/src/config.rs b/crates/within-py/src/config.rs index 7e280ba2..8b3068cf 100644 --- a/crates/within-py/src/config.rs +++ b/crates/within-py/src/config.rs @@ -55,8 +55,6 @@ impl PyApproxCholConfig { #[pyclass(frozen, module = "within._within")] #[pyo3(name = "ApproxSchurConfig")] pub struct PyApproxSchurConfig { - #[pyo3(get)] - pub seed: u64, #[pyo3(get)] pub split: u32, } @@ -64,23 +62,20 @@ pub struct PyApproxSchurConfig { #[pymethods] impl PyApproxSchurConfig { #[new] - #[pyo3(signature = (seed=0, split=1))] - fn new(seed: u64, split: u32) -> PyResult { + #[pyo3(signature = (split=1))] + fn new(split: u32) -> PyResult { if split == 0 { return Err(pyo3::exceptions::PyValueError::new_err( "split must be >= 1", )); } - Ok(Self { seed, split }) + Ok(Self { split }) } } impl PyApproxSchurConfig { pub(crate) fn to_native(&self) -> ApproxSchurConfig { - ApproxSchurConfig { - seed: self.seed, - split: self.split, - } + ApproxSchurConfig { split: self.split } } } @@ -117,7 +112,7 @@ impl PySchur { fn __repr__(&self) -> String { match &self.inner { SchurMode::Approximate(cfg) => { - format!("Schur.approximate(seed={}, split={})", cfg.seed, cfg.split) + format!("Schur.approximate(split={})", cfg.split) } SchurMode::Exact => "Schur.exact()".to_string(), } diff --git a/crates/within/Cargo.toml b/crates/within/Cargo.toml index df6bb103..a4490c16 100644 --- a/crates/within/Cargo.toml +++ b/crates/within/Cargo.toml @@ -12,7 +12,7 @@ categories.workspace = true [dependencies] schwarz-precond = { workspace = true, features = ["serde"] } -approx-chol = { version = "0.6.0", features = ["serde"] } +approx-chol = { version = "0.7.0", features = ["serde"] } serde = { workspace = true, features = ["rc"] } postcard.workspace = true ndarray = "0.17" diff --git a/crates/within/src/block_elim.rs b/crates/within/src/block_elim.rs deleted file mode 100644 index ce901553..00000000 --- a/crates/within/src/block_elim.rs +++ /dev/null @@ -1,17 +0,0 @@ -//! Bipartite block-elimination local solver for Schwarz factor-pair subdomains. -//! -//! [`BlockElimSolver`] handles a 2x2 bipartite Gramian `[D_q, C; C^T, D_r]` by -//! eliminating the larger diagonal block, then solving the Schur complement -//! via approximate or dense Cholesky (see [`schur`] and [`factor`]). - -pub(crate) mod csr_matrix; -pub(crate) mod factor; -pub(crate) mod schur; -pub(crate) mod solver; - -pub(crate) use solver::BlockElimSolver; - -pub(crate) use crate::linalg::compensated_sum; - -#[cfg(test)] -mod tests; diff --git a/crates/within/src/block_elim/csr_matrix.rs b/crates/within/src/block_elim/csr_matrix.rs deleted file mode 100644 index 36843613..00000000 --- a/crates/within/src/block_elim/csr_matrix.rs +++ /dev/null @@ -1,82 +0,0 @@ -//! Validated square `f64`/`u32` matrix in Compressed Sparse Row (CSR) format. -//! -//! Holds the reduced Schur-complement / Laplacian system that -//! [`schur`](super::schur) assembles and [`factor`](super::factor) hands to the -//! sparse factor backend — a `block_elim`-internal representation. It is -//! distinct from the crate-level [`CsrBlock`](crate::csr_block::CsrBlock), which -//! carries the rectangular cross-factor coupling block produced in `domain` and -//! consumed here. - -/// Block elimination produces the CSR invariants, so [`CsrMatrix::new`] only debug-checks them. -#[derive(Clone)] -pub(crate) struct CsrMatrix { - indptr: Vec, - indices: Vec, - data: Vec, - n: usize, -} - -impl CsrMatrix { - /// Trusted in release, `assert!`-checked in debug. - pub(crate) fn new(indptr: Vec, indices: Vec, data: Vec, n: usize) -> Self { - #[cfg(debug_assertions)] - { - assert_eq!(indptr.len(), n + 1, "indptr length must equal n + 1"); - assert_eq!(indptr[0], 0, "indptr[0] must be 0"); - assert_eq!( - indices.len(), - data.len(), - "indices and data must have equal length" - ); - assert_eq!( - indices.len(), - indptr[n] as usize, - "indices/data length must equal indptr[n]" - ); - for row in 0..n { - let start = indptr[row] as usize; - let end = indptr[row + 1] as usize; - assert!(end >= start, "indptr must be non-decreasing"); - let cols = &indices[start..end]; - assert!( - cols.iter().all(|&col| (col as usize) < n), - "column index out of bounds" - ); - assert!( - cols.windows(2).all(|w| w[0] < w[1]), - "row column indices must be strictly ascending" - ); - } - } - Self { - indptr, - indices, - data, - n, - } - } - - /// Row pointer array (length `n + 1`). - #[inline] - pub(crate) fn indptr(&self) -> &[u32] { - &self.indptr - } - - /// Column indices (length `nnz`). - #[inline] - pub(crate) fn indices(&self) -> &[u32] { - &self.indices - } - - /// Non-zero values (length `nnz`), parallel to `indices`. - #[inline] - pub(crate) fn data(&self) -> &[f64] { - &self.data - } - - /// Matrix dimension (the matrix is always square). - #[inline] - pub(crate) fn n(&self) -> usize { - self.n - } -} diff --git a/crates/within/src/block_elim/factor.rs b/crates/within/src/block_elim/factor.rs deleted file mode 100644 index 76d5c6b0..00000000 --- a/crates/within/src/block_elim/factor.rs +++ /dev/null @@ -1,248 +0,0 @@ -//! Reduced-system factor for Schur-complement local solves. -//! -//! [`ReducedFactor`] wraps an `approx-chol` factor, either directly or behind a -//! Gremban cover. [`factor_sparse`] bridges to the `approx-chol` builder. - -use approx_chol::low_level::Builder; -use approx_chol::{CsrRef, Factor}; -use schwarz_precond::LocalSolveError; - -use super::csr_matrix::CsrMatrix; -use crate::domain::Grounding; -use crate::BuildError; - -/// Reduced-system factor for Schur-complement local solves. -#[derive(Clone, serde::Serialize)] -pub(crate) enum ReducedFactor { - // Postcard encodes discriminants by declaration order and the fixture pins Direct = 0. - /// Factor of the reduced Schur CSR itself, carrying the gauge applied around the solve. - Direct { - /// Factor of the reduced Schur. - factor: Factor, - /// Gauge of the reduced system. - grounding: Grounding, - }, - /// Gremban cover of a signed reduced Schur, kept here so the operator stays single-sized. - Cover { - /// Factor of the doubled cover; the tail past `2*m` is the grounding augmentation. - inner: Factor, - /// Single signed reduced dimension exposed to the caller. - m: usize, - }, -} - -impl ReducedFactor { - /// The gauge the operator-level solve applies; `None` for a self-grounding cover. - pub(crate) fn grounding(&self) -> Option { - match self { - Self::Direct { grounding, .. } => Some(*grounding), - Self::Cover { .. } => None, - } - } - - /// Explicit ground vertex appended past the kept block; `None` for the kept block alone. - pub(crate) fn explicit_ground_index(&self, n_kept: usize) -> Option { - (self.grounding() == Some(Grounding::Grounded) && self.input_dimension() == n_kept + 1) - .then_some(n_kept) - } - - /// Whether the factor was built over a kept block of this size, with or without ground. - pub(crate) fn spans_kept_block(&self, n_kept: usize) -> bool { - self.input_dimension() == n_kept || self.explicit_ground_index(n_kept).is_some() - } - - /// Dimension handed to the backend, which adds grounding vertices of its own. - pub(crate) fn input_dimension(&self) -> usize { - match self { - Self::Direct { factor, .. } => factor.original_n(), - // The cover is hidden behind the single signed interface. - Self::Cover { m, .. } => *m, - } - } - - pub(crate) fn solve_dimension(&self) -> usize { - match self { - Self::Direct { factor, .. } => factor.n(), - Self::Cover { m, .. } => *m, - } - } - - /// Length of the larger system only [`Self::Cover`] embeds into. - pub(crate) fn embed_len(&self) -> usize { - match self { - Self::Direct { .. } => 0, - Self::Cover { inner, .. } => inner.n(), - } - } - - /// `x` spans [`Self::solve_dimension`]; `embed` is at least [`Self::embed_len`] long. - pub(crate) fn solve_in_place( - &self, - x: &mut [f64], - embed: &mut [f64], - ) -> Result<(), LocalSolveError> { - match self { - Self::Direct { factor, .. } => { - debug_assert_eq!(factor.n(), x.len()); - solve_approx(factor, x) - } - Self::Cover { inner, m } => { - debug_assert_eq!(*m, x.len()); - let cover_n = inner.n(); - let buf = &mut embed[..cover_n]; - // Grounding vertices past the `2m` cover nodes have zero RHS, so clear the tail. - buf[..*m].copy_from_slice(x); - for (out, &v) in buf[*m..2 * *m].iter_mut().zip(x.iter()) { - *out = -v; - } - for slot in buf[2 * *m..].iter_mut() { - *slot = 0.0; - } - solve_approx(inner, buf)?; - // Read back the antisymmetric solution: x = (x⁺ - x⁻) / 2. - for (i, out) in x.iter_mut().enumerate() { - *out = 0.5 * (buf[i] - buf[*m + i]); - } - Ok(()) - } - } - } -} - -/// Wire mirror of [`ReducedFactor`]; a bare inner [`Factor`] makes `Cover`-of-`Cover` undecodable. -#[derive(serde::Deserialize)] -enum ReducedFactorWire { - Direct { - factor: Factor, - grounding: Grounding, - }, - Cover { - inner: Factor, - m: usize, - }, -} - -impl<'de> serde::Deserialize<'de> for ReducedFactor { - fn deserialize>(deserializer: D) -> Result { - use serde::de::Error; - match ReducedFactorWire::deserialize(deserializer)? { - ReducedFactorWire::Direct { factor, grounding } => { - Ok(ReducedFactor::Direct { factor, grounding }) - } - ReducedFactorWire::Cover { inner, m } => { - // Outside `[2m, 2m + 2]` the antisymmetric embed over- or under-runs. - let two_m = m - .checked_mul(2) - .ok_or_else(|| D::Error::custom("Cover dimension m too large"))?; - let cover_dim = inner.n(); - if cover_dim < two_m || cover_dim > two_m.saturating_add(2) { - return Err(D::Error::custom( - "Cover inner factor dimension inconsistent with m", - )); - } - Ok(ReducedFactor::Cover { inner, m }) - } - } - } -} - -fn solve_approx(f: &Factor, x: &mut [f64]) -> Result<(), LocalSolveError> { - f.solve_in_place(x) - .map_err(|e| LocalSolveError::BackendFailed { - context: "within.local.block_elim.reduced_approx", - message: e.to_string(), - }) -} - -/// Returns the `approx-chol` error unmapped, so the caller can spot an unusable exact pivot. -pub(crate) fn factor_sparse( - matrix: &CsrMatrix, - config: approx_chol::Config, -) -> Result { - let csr = CsrRef::new( - matrix.indptr(), - matrix.indices(), - matrix.data(), - u32::try_from(matrix.n()).expect("Schur complement dimension exceeds u32::MAX"), - )?; - Builder::new(config).build(csr) -} - -pub(crate) fn local_solver_build(e: approx_chol::Error) -> BuildError { - BuildError::LocalSolverBuild(format!("failed Schur complement factorization: {e}")) -} - -#[cfg(test)] -mod tests { - use super::*; - use crate::config::{ApproxCholConfig, DEFAULT_DENSE_SCHUR_THRESHOLD}; - use approx_chol::ExactFailure; - - #[test] - fn cover_reduced_factor_solves_signed_system() { - // Pins the `ReducedFactor::Cover` embed/read-back, not approx-chol's accuracy. - let cover = CsrMatrix::new( - vec![0, 2, 4, 6, 8], - vec![0, 3, 1, 2, 1, 2, 0, 3], - vec![2.0, -1.0, 2.0, -1.0, -1.0, 2.0, -1.0, 2.0], - 4, - ); - let config = ApproxCholConfig { - seed: 0, - split_merge: Some(2), - } - .to_approx_chol(DEFAULT_DENSE_SCHUR_THRESHOLD, ExactFailure::Error); - let inner = factor_sparse(&cover, config).expect("factor cover"); - let reduced = ReducedFactor::Cover { inner, m: 2 }; - assert_eq!(reduced.input_dimension(), 2); - assert_eq!(reduced.solve_dimension(), 2); - - let b = [1.0, 0.5]; - let mut x = b; - let mut embed = vec![0.0; reduced.embed_len()]; - reduced - .solve_in_place(&mut x, &mut embed) - .expect("cover solve"); - - // Residual of the signed system M x = b (M = [[2,1],[1,2]]). - let r0 = 2.0 * x[0] + x[1] - b[0]; - let r1 = x[0] + 2.0 * x[1] - b[1]; - assert!(r0.hypot(r1) < 1e-9, "residual too large: ({r0}, {r1})"); - } - - /// A 2×2 SDDM system; approx-chol grounds its surplus, so the factor is - /// dimension 3 over an input of 2. - fn approx_2x2() -> Factor { - let m = CsrMatrix::new( - vec![0, 2, 4], - vec![0, 1, 0, 1], - vec![2.0, -1.0, -1.0, 2.0], - 2, - ); - let config = ApproxCholConfig::default() - .to_approx_chol(DEFAULT_DENSE_SCHUR_THRESHOLD, ExactFailure::Error); - factor_sparse(&m, config).expect("factor 2x2") - } - - #[test] - fn valid_direct_round_trips() { - let bytes = postcard::to_stdvec(&ReducedFactor::Direct { - factor: approx_2x2(), - grounding: Grounding::Floating, - }) - .expect("serialize"); - let restored: ReducedFactor = postcard::from_bytes(&bytes).expect("deserialize"); - assert_eq!(restored.input_dimension(), 2); - } - - #[test] - fn cover_with_undersized_inner_is_rejected() { - // The inner factor cannot hold the antisymmetric embed for m = 5, which needs 10 nodes. - let bad = ReducedFactor::Cover { - inner: approx_2x2(), - m: 5, - }; - let bytes = postcard::to_stdvec(&bad).expect("serialize"); - assert!(postcard::from_bytes::(&bytes).is_err()); - } -} diff --git a/crates/within/src/block_elim/schur.rs b/crates/within/src/block_elim/schur.rs deleted file mode 100644 index 8b64ba93..00000000 --- a/crates/within/src/block_elim/schur.rs +++ /dev/null @@ -1,473 +0,0 @@ -//! Schur complement of the eliminated diagonal block in a bipartite SDDM -//! system — SDDM in, SDDM on the kept levels out. -//! -//! [`exact`] accumulates the true reduced rows; [`sampled`] approximates each -//! eliminated star's clique on the explicit augmented graph. - -use super::compensated_sum; -use super::csr_matrix::CsrMatrix; -use approx_chol::low_level::CliqueTreeSampler; -use rayon::prelude::*; - -use crate::config::ApproxSchurConfig; -use crate::csr_block::to_u32; -use crate::domain::{Grounding, SddmMatrix}; - -/// Whether the kept rows are split across threads. The subdomain loop reaching here is -/// already parallel, so splitting pays only once the reduction itself is large enough to -/// outweigh the tasks it adds for that loop to steal. -#[derive(Debug, Clone, Copy, PartialEq, Eq)] -pub(crate) enum RowSplit { - Sequential, - Parallel, -} - -struct UpperSchurRow { - columns: Vec, - values: Vec, -} - -/// Multiply-adds below which splitting the kept rows costs more than it saves. -const PAR_ROW_SPLIT_THRESHOLD: usize = 100_000; - -/// Exact Schur complement, accumulated per keep-row without materializing intermediate edges. -pub(crate) fn exact( - matrix: &SddmMatrix, - inv_diagonal_eliminated: &[f64], - split: RowSplit, -) -> CsrMatrix { - let n_keep = matrix.n_kept(); - // `extract_sparse_row` zeroes what it read, so one workspace serves a whole run. - let row = |i: usize, work: &mut Vec, touched: &mut Vec| { - compute_schur_row_dense(matrix, inv_diagonal_eliminated, i, work, touched); - let result = extract_sparse_row(i, work, touched); - touched.clear(); - result - }; - - let rows: Vec = match split { - RowSplit::Parallel => (0..n_keep) - .into_par_iter() - .map_init( - || (vec![0.0f64; n_keep], Vec::new()), - |(work, touched), i| row(i, work, touched), - ) - .collect(), - RowSplit::Sequential => { - let mut work = vec![0.0f64; n_keep]; - let mut touched = Vec::new(); - (0..n_keep) - .map(|i| row(i, &mut work, &mut touched)) - .collect() - } - }; - - assemble_schur_csr(rows, n_keep) -} - -/// Sampled Schur complement as a Laplacian; a ground vertex stays an ordinary final vertex. -pub(crate) fn sampled(matrix: &SddmMatrix, config: &ApproxSchurConfig) -> CsrMatrix { - let edges = par_emit(matrix, config); - let n = matrix.n_kept() + usize::from(matrix.grounding == Grounding::Grounded); - build_laplacian_csr(&edges, n) -} - -pub(crate) fn exact_for_factor(matrix: &SddmMatrix, inv_diagonal_eliminated: &[f64]) -> CsrMatrix { - let split = if exact_flops(matrix) < PAR_ROW_SPLIT_THRESHOLD { - RowSplit::Sequential - } else { - RowSplit::Parallel - }; - let principal = exact(matrix, inv_diagonal_eliminated, split); - let surplus = reduced_surplus(matrix, inv_diagonal_eliminated); - build_explicit_laplacian(&principal, &surplus, matrix.grounding) -} - -/// Upper-triangle multiply-adds: `Σ_k nnz(k)(nnz(k)+1)/2`. -fn exact_flops(matrix: &SddmMatrix) -> usize { - matrix - .cross_tab - .c - .indptr - .windows(2) - .map(|bounds| { - let width = (bounds[1] - bounds[0]) as usize; - width * (width + 1) / 2 - }) - .sum() -} - -/// Scatter the diagonal-and-upper part of Schur row `i`, recording touched columns. -fn compute_schur_row_dense( - matrix: &SddmMatrix, - inv_diagonal_eliminated: &[f64], - i: usize, - work: &mut [f64], - touched: &mut Vec, -) { - work[i] = matrix.diagonal_kept()[i]; - touched.push(i); - - let elim_to_keep = &matrix.cross_tab.c; - for (k, w) in matrix.cross_tab.ct().row(i) { - let inv = inv_diagonal_eliminated[k]; - let start = elim_to_keep.indptr[k] as usize; - let end = elim_to_keep.indptr[k + 1] as usize; - let columns = &elim_to_keep.indices[start..end]; - let upper_start = columns.partition_point(|&j| (j as usize) < i); - for position in start + upper_start..end { - let j = elim_to_keep.indices[position] as usize; - let v = elim_to_keep.data[position]; - if work[j] == 0.0 && j != i { - touched.push(j); - } - work[j] -= w * v * inv; - } - } -} - -/// Extract non-zeros, keeping a numerically-zero diagonal so SDDM structure survives. -fn extract_sparse_row(i: usize, work: &mut [f64], touched: &mut [usize]) -> UpperSchurRow { - touched.sort_unstable(); - let mut columns = Vec::with_capacity(touched.len()); - let mut values = Vec::with_capacity(touched.len()); - for &j in touched.iter() { - let v = work[j]; - if v != 0.0 || j == i { - columns.push(to_u32(j)); - values.push(v); - } - work[j] = 0.0; - } - UpperSchurRow { columns, values } -} - -/// Assemble full rows from diagonal-and-upper parts; the strict upper mirrors into the lower. -fn assemble_schur_csr(rows: Vec, n_keep: usize) -> CsrMatrix { - let mut row_offsets = vec![0u32; n_keep + 1]; - for (i, row) in rows.iter().enumerate() { - row_offsets[i + 1] += to_u32(row.columns.len()); - for &j in &row.columns { - if j as usize != i { - row_offsets[j as usize + 1] += 1; - } - } - } - for i in 0..n_keep { - row_offsets[i + 1] += row_offsets[i]; - } - let total_nnz = row_offsets[n_keep] as usize; - let mut column_indices = vec![0u32; total_nnz]; - let mut values = vec![0.0f64; total_nnz]; - - // Row layout: [mirrored lower | diagonal and upper]; ascending source rows keep it sorted. - let mut lower_write_positions = row_offsets[..n_keep].to_vec(); - for (i, row) in rows.into_iter().enumerate() { - let upper_start = row_offsets[i + 1] as usize - row.columns.len(); - column_indices[upper_start..upper_start + row.columns.len()].copy_from_slice(&row.columns); - values[upper_start..upper_start + row.values.len()].copy_from_slice(&row.values); - for (&j, &v) in row.columns.iter().zip(&row.values) { - if j as usize != i { - let position = lower_write_positions[j as usize] as usize; - column_indices[position] = to_u32(i); - values[position] = v; - lower_write_positions[j as usize] += 1; - } - } - } - CsrMatrix::new(row_offsets, column_indices, values, n_keep) -} - -/// Edges sorted by `(lo, hi)` land both triangles in column order without per-row sorting. -fn build_laplacian_csr(edges: &[Edge], n: usize) -> CsrMatrix { - debug_assert!(edges.iter().all(|&(lo, hi, _)| lo < hi)); - - let mut lower_count = vec![0u32; n]; - let mut upper_count = vec![0u32; n]; - let mut diag = vec![0.0; n]; - for &(lo, hi, w) in edges { - diag[lo as usize] += w; - upper_count[lo as usize] += 1; // row lo gets col hi (upper) - lower_count[hi as usize] += 1; // row hi gets col lo (lower) - diag[hi as usize] += w; - } - - // Row layout: [lower entries | diagonal | upper entries]. - let mut offsets = vec![0u32; n + 1]; - for i in 0..n { - offsets[i + 1] = offsets[i] + lower_count[i] + 1 + upper_count[i]; - } - let total_nnz = offsets[n] as usize; - let mut indices = vec![0u32; total_nnz]; - let mut data = vec![0.0f64; total_nnz]; - - let mut lower_cursor: Vec = (0..n).map(|i| offsets[i]).collect(); - let mut upper_cursor: Vec = (0..n).map(|i| offsets[i] + lower_count[i] + 1).collect(); - for i in 0..n { - let pos = (offsets[i] + lower_count[i]) as usize; - indices[pos] = to_u32(i); - data[pos] = diag[i]; - } - - for &(lo, hi, w) in edges { - let lo_idx = lo as usize; - let hi_idx = hi as usize; - // Upper triangle: row lo, column hi. - let pos = upper_cursor[lo_idx] as usize; - indices[pos] = hi; - data[pos] = -w; - upper_cursor[lo_idx] += 1; - // Lower triangle: row hi, column lo. - let pos = lower_cursor[hi_idx] as usize; - indices[pos] = lo; - data[pos] = -w; - lower_cursor[hi_idx] += 1; - } - - CsrMatrix::new(offsets, indices, data, n) -} - -fn reduced_surplus(matrix: &SddmMatrix, inv_diagonal_eliminated: &[f64]) -> Vec { - if matrix.grounding == Grounding::Floating { - return vec![0.0; matrix.n_kept()]; - } - let scaled: Vec = inv_diagonal_eliminated - .iter() - .zip(matrix.surplus_eliminated()) - .map(|(&inv_diag, &surplus)| inv_diag * surplus) - .collect(); - let mut surplus = vec![0.0; matrix.n_kept()]; - matrix - .cross_tab - .ct() - .spmv_assign_add(&scaled, matrix.surplus_kept(), &mut surplus, false); - surplus -} - -pub(super) fn build_explicit_laplacian( - principal: &CsrMatrix, - surplus: &[f64], - grounding: Grounding, -) -> CsrMatrix { - let n_keep = principal.n(); - let ground = to_u32(n_keep); - let grounded = grounding == Grounding::Grounded; - let n = n_keep + usize::from(grounded); - let mut indptr = Vec::with_capacity(n + 1); - let mut indices = Vec::new(); - let mut data = Vec::new(); - indptr.push(0); - - for (i, &row_surplus) in surplus.iter().enumerate().take(n_keep) { - let start = principal.indptr()[i] as usize; - let end = principal.indptr()[i + 1] as usize; - let mut adjacency = 0.0; - let mut diagonal_position = None; - for (&j, &value) in principal.indices()[start..end] - .iter() - .zip(&principal.data()[start..end]) - { - if j as usize == i { - diagonal_position = Some(indices.len()); - indices.push(j); - data.push(0.0); - } else { - adjacency -= value; - indices.push(j); - data.push(value); - } - } - let diagonal_position = diagonal_position.expect("exact Schur row must contain a diagonal"); - data[diagonal_position] = adjacency + row_surplus; - if grounded && row_surplus > 0.0 { - indices.push(ground); - data.push(-row_surplus); - } - indptr.push(to_u32(indices.len())); - } - - if grounded { - for (i, &value) in surplus.iter().enumerate() { - if value > 0.0 { - indices.push(to_u32(i)); - data.push(-value); - } - } - indices.push(ground); - data.push(compensated_sum(surplus)); - indptr.push(to_u32(indices.len())); - } - CsrMatrix::new(indptr, indices, data, n) -} - -/// Undirected fill edge: `(lo_col, hi_col, weight)` with `lo_col < hi_col`. -type Edge = (u32, u32, f64); - -/// One eliminated vertex's keep-block neighbors, referencing the cross-tab CSR for zero copy. -struct Star<'a> { - /// Eliminated vertex index (used for deterministic seeding). - index: usize, - /// Neighbor columns in the keep-block. - col_indices: &'a [u32], - /// Edge weights to each neighbor. - weights: &'a [f64], -} - -impl Star<'_> { - fn degree(&self) -> usize { - self.col_indices.len() - } -} - -/// One rayon task's sampling state. The sampler carries scratch of its own, so it is -/// built per task alongside the buffers rather than per star. -struct StarWorkspace { - edges: Vec, - entries: Vec<(u32, f64)>, - sampler: CliqueTreeSampler, -} - -impl StarWorkspace { - fn new(config: &ApproxSchurConfig) -> Self { - Self { - edges: Vec::new(), - entries: Vec::new(), - sampler: CliqueTreeSampler::new(config.seed, Some(config.split)), - } - } - - /// Ground surplus is one more incident edge, so a sampled star's capacity is its diagonal. - fn sample_star(&mut self, star: &Star, ground: Option<(u32, f64)>) { - self.entries.clear(); - for (&col, &w) in star.col_indices.iter().zip(star.weights) { - self.entries.push((col, w)); - } - if let Some((ground, surplus)) = ground { - if surplus > 0.0 { - self.entries.push((ground, surplus)); - } - } - if self.entries.len() <= 1 { - return; - } - self.sampler - .sample(star.index as u64, &self.entries, &mut self.edges); - } -} - -/// Create a zero-copy [`Star`] view for eliminated vertex `k`. -fn star(matrix: &SddmMatrix, k: usize) -> Star<'_> { - let elim_to_keep = &matrix.cross_tab.c; - let start = elim_to_keep.indptr[k] as usize; - let end = elim_to_keep.indptr[k + 1] as usize; - Star { - index: k, - col_indices: &elim_to_keep.indices[start..end], - weights: &elim_to_keep.data[start..end], - } -} - -fn par_emit(matrix: &SddmMatrix, config: &ApproxSchurConfig) -> Vec { - // The total order fixes per-`(lo, hi)` summation, so the result is scheduling-independent. - let n_kept = matrix.n_kept(); - let surplus_eliminated = matrix.surplus_eliminated(); - let ground_vertex = (matrix.grounding == Grounding::Grounded) - .then(|| u32::try_from(n_kept).expect("ground vertex exceeds u32::MAX")); - let mut chunks = (0..matrix.n_eliminated()) - .into_par_iter() - .fold( - || StarWorkspace::new(config), - |mut work, k| { - let star = star(matrix, k); - if star.degree() > 0 { - let ground = ground_vertex.map(|g| (g, surplus_eliminated[k])); - work.sample_star(&star, ground); - } - work - }, - ) - .map(|work| work.edges) - .collect::>(); - if let Some(ground) = ground_vertex { - chunks.push( - matrix - .surplus_kept() - .iter() - .enumerate() - .filter(|(_, surplus)| **surplus > 0.0) - .map(|(i, &surplus)| (i as u32, ground, surplus)) - .collect(), - ); - } - sort_and_dedup(chunks, n_kept) -} - -/// Sort into total `(lo, hi, weight)` order; the weight tiebreak makes the Schur reproducible. -fn sort_and_dedup(chunks: Vec>, n_kept: usize) -> Vec { - let total = chunks.iter().map(Vec::len).sum::(); - if total <= 1 { - return chunks.concat(); - } - // Same order either way; under two edges per kept vertex the parallel comparison sort wins. - let dense = total >= 2 * n_kept; - let mut edges = if dense { - counting_sort_by_lo(&chunks, n_kept) - } else { - chunks.concat() - }; - drop(chunks); - if dense { - let by_hi_weight = |a: &Edge, b: &Edge| a.1.cmp(&b.1).then_with(|| a.2.total_cmp(&b.2)); - edges - .par_chunk_by_mut(|a, b| a.0 == b.0) - .for_each(|run| run.sort_unstable_by(by_hi_weight)); - } else { - edges.par_sort_unstable_by(|a, b| { - a.0.cmp(&b.0) - .then(a.1.cmp(&b.1)) - .then_with(|| a.2.total_cmp(&b.2)) - }); - } - let mut write = 0; - for read in 1..edges.len() { - if edges[write].0 == edges[read].0 && edges[write].1 == edges[read].1 { - edges[write].2 += edges[read].2; - } else { - write += 1; - edges[write] = edges[read]; - } - } - edges.truncate(write + 1); - edges -} - -/// Stable counting sort by `lo` in O(E + n_kept), scattered once from the per-task chunks. -fn counting_sort_by_lo(chunks: &[Vec], n_kept: usize) -> Vec { - let mut cursors = vec![0usize; n_kept + 1]; - // Nested loops, not `flatten()`: the flat iterator cost ~10% on both passes at 1M+ edges. - for chunk in chunks { - for e in chunk { - // `lo < n_kept` always holds: the ground vertex has the maximum id and lands on `hi`. - debug_assert!( - (e.0 as usize) < n_kept, - "counting sort key `lo` must be a kept-block id (< n_kept)" - ); - cursors[e.0 as usize + 1] += 1; - } - } - for i in 1..cursors.len() { - cursors[i] += cursors[i - 1]; - } - let mut out = vec![(0u32, 0u32, 0.0f64); cursors[n_kept]]; - for chunk in chunks { - for &e in chunk { - let cursor = &mut cursors[e.0 as usize]; - out[*cursor] = e; - *cursor += 1; - } - } - out -} - -#[cfg(test)] -mod tests; diff --git a/crates/within/src/block_elim/schur/tests.rs b/crates/within/src/block_elim/schur/tests.rs deleted file mode 100644 index 33f00701..00000000 --- a/crates/within/src/block_elim/schur/tests.rs +++ /dev/null @@ -1,205 +0,0 @@ -use rstest::rstest; - -use super::*; - -fn sparse_to_dense(matrix: &CsrMatrix) -> Vec> { - let n = matrix.n(); - let mut dense = vec![vec![0.0; n]; n]; - for (i, row) in dense.iter_mut().enumerate().take(n) { - let start = matrix.indptr()[i] as usize; - let end = matrix.indptr()[i + 1] as usize; - for idx in start..end { - let j = matrix.indices()[idx] as usize; - row[j] = matrix.data()[idx]; - } - } - dense -} - -fn dense_exact_schur( - c_dense: &[f64], - n_rows: usize, - n_cols: usize, - row_diag: &[f64], - col_diag: &[f64], - eliminate_rows: bool, -) -> Vec> { - if eliminate_rows { - let mut s = vec![vec![0.0; n_cols]; n_cols]; - for i in 0..n_cols { - s[i][i] = col_diag[i]; - } - for k in 0..n_rows { - let inv = if row_diag[k] > 0.0 { - 1.0 / row_diag[k] - } else { - 0.0 - }; - for i in 0..n_cols { - let cki = c_dense[k * n_cols + i]; - for j in 0..n_cols { - let ckj = c_dense[k * n_cols + j]; - s[i][j] -= cki * inv * ckj; - } - } - } - s - } else { - let mut s = vec![vec![0.0; n_rows]; n_rows]; - for i in 0..n_rows { - s[i][i] = row_diag[i]; - } - for k in 0..n_cols { - let inv = if col_diag[k] > 0.0 { - 1.0 / col_diag[k] - } else { - 0.0 - }; - for i in 0..n_rows { - let cik = c_dense[i * n_cols + k]; - for j in 0..n_rows { - let cjk = c_dense[j * n_cols + k]; - s[i][j] -= cik * inv * cjk; - } - } - } - s - } -} - -fn assert_dense_close(lhs: &[Vec], rhs: &[Vec], tol: f64) { - assert_eq!(lhs.len(), rhs.len(), "row count mismatch"); - for i in 0..lhs.len() { - assert_eq!(lhs[i].len(), rhs[i].len(), "col count mismatch on row {i}"); - for j in 0..lhs[i].len() { - assert!( - (lhs[i][j] - rhs[i][j]).abs() <= tol, - "mismatch at ({i}, {j}): lhs={}, rhs={}", - lhs[i][j], - rhs[i][j] - ); - } - } -} - -#[test] -fn exact_schur_matches_dense_reference_when_eliminating_rows() { - // C is 3x2, so the row block is eliminated (n_rows >= n_cols). - let c_dense = vec![1.0, 2.0, 3.0, 0.0, 0.0, 4.0]; - let diagonal = vec![5.0, 6.0, 8.0, 7.0, 9.0]; - let (row_diag, col_diag) = diagonal.split_at(3); - let inv_diagonal: Vec = row_diag.iter().map(|d| 1.0 / d).collect(); - let expected = dense_exact_schur(&c_dense, 3, 2, row_diag, col_diag, true); - - let matrix = SddmMatrix::from_dense_for_test(&c_dense, 3, 2, diagonal, Grounding::Grounded); - let got = sparse_to_dense(&exact(&matrix, &inv_diagonal, RowSplit::Parallel)); - assert_dense_close(&got, &expected, 1e-12); -} - -/// The sequential arm reuses one scatter workspace across rows, so a row that failed to -/// reset it would read the previous row's entries. -#[rstest] -fn both_row_splits_produce_the_same_complement( - #[values(Grounding::Grounded, Grounding::Floating)] grounding: Grounding, -) { - let c_dense = vec![1.0, 2.0, 3.0, 0.0, 0.0, 4.0]; - let diagonal = vec![5.0, 6.0, 8.0, 7.0, 9.0]; - let (row_diag, _) = diagonal.split_at(3); - let inv_diagonal: Vec = row_diag.iter().map(|d| 1.0 / d).collect(); - - let operator = SddmMatrix::from_dense_for_test(&c_dense, 3, 2, diagonal, grounding); - let parallel = exact(&operator, &inv_diagonal, RowSplit::Parallel); - let sequential = exact(&operator, &inv_diagonal, RowSplit::Sequential); - - assert_eq!(sequential.indptr(), parallel.indptr()); - assert_eq!(sequential.indices(), parallel.indices()); - assert_eq!(sequential.data(), parallel.data()); -} - -#[test] -fn approximate_schur_is_seed_deterministic_and_laplacian_like() { - // Stating the sums the traversal must reproduce keeps the Laplacian assertion an independent check. - let c_dense = vec![1.0, 2.0, 3.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0]; - let matrix = SddmMatrix::from_dense_for_test( - &c_dense, - 3, - 3, - vec![6.0, 1.0, 1.0, 2.0, 3.0, 3.0], - Grounding::Floating, - ); - let config = crate::config::ApproxSchurConfig { - seed: 12345, - ..Default::default() - }; - - let a = sampled(&matrix, &config); - let b = sampled(&matrix, &config); - - assert_eq!(a.indptr(), b.indptr()); - assert_eq!(a.indices(), b.indices()); - assert_eq!(a.data(), b.data()); - - let dense = sparse_to_dense(&a); - for (i, row) in dense.iter().enumerate() { - let mut row_sum = 0.0; - for (j, &value) in row.iter().enumerate() { - row_sum += value; - assert!( - (value - dense[j][i]).abs() <= 1e-12, - "matrix not symmetric at ({i}, {j})" - ); - if i != j { - assert!(value <= 1e-12, "off-diagonal should be non-positive"); - } - } - assert!(row_sum.abs() <= 1e-10, "row {i} sum is not near zero"); - assert!(row[i] >= -1e-12, "diagonal should be non-negative"); - } -} - -#[test] -fn sampled_schur_carries_surplus_exactly_on_low_degree_stars() { - // Every eliminated star has at most two entries including ground, so sampling is exact. - let c_dense = vec![1.0, 0.0, 2.0, 3.0, 0.0, 4.0]; - let diagonal = vec![3.0, 5.0, 6.0, 3.5, 7.0]; - let (row_diag, col_diag) = diagonal.split_at(3); - let expected = dense_exact_schur(&c_dense, 3, 2, row_diag, col_diag, true); - let matrix = SddmMatrix::from_dense_for_test(&c_dense, 3, 2, diagonal, Grounding::Grounded); - - let sampled_dense = sparse_to_dense(&sampled(&matrix, &Default::default())); - let sampled_principal: Vec> = sampled_dense[..expected.len()] - .iter() - .map(|row| row[..expected.len()].to_vec()) - .collect(); - assert_dense_close(&sampled_principal, &expected, 1e-12); - - // Surplus becomes an explicit ground vertex, so the augmented matrix is a Laplacian. - for (i, row) in sampled_dense.iter().enumerate() { - let row_sum: f64 = row.iter().sum(); - assert!(row_sum.abs() < 1e-12, "row {i} sum is {row_sum}"); - } -} - -#[test] -fn assemble_schur_csr_mirrors_upper_rows() { - let rows = vec![ - UpperSchurRow { - columns: vec![0, 2], - values: vec![4.0, -1.25], - }, - UpperSchurRow { - columns: vec![1, 2], - values: vec![5.0, -2.5], - }, - UpperSchurRow { - columns: vec![2], - values: vec![6.0], - }, - ]; - - let matrix = assemble_schur_csr(rows, 3); - - assert_eq!(matrix.indptr(), &[0, 2, 4, 7]); - assert_eq!(matrix.indices(), &[0, 2, 1, 2, 0, 1, 2]); - assert_eq!(matrix.data(), &[4.0, -1.25, 5.0, -2.5, -1.25, -2.5, 6.0]); -} diff --git a/crates/within/src/block_elim/solver.rs b/crates/within/src/block_elim/solver.rs deleted file mode 100644 index 39924dd4..00000000 --- a/crates/within/src/block_elim/solver.rs +++ /dev/null @@ -1,468 +0,0 @@ -use std::borrow::Borrow; -use std::sync::Arc; - -use approx_chol::{ExactFailure, Factor}; -use rayon::prelude::*; -use schwarz_precond::{LocalSolveError, LocalSolver}; - -use crate::config::{LocalSolverConfig, SchurMode}; -use crate::csr_block::{CsrBlock, PAR_SPMV_THRESHOLD}; -use crate::domain::{CoordinateMap, CrossTab, Grounding, LocalComponent, MatrixForm, SddmMatrix}; -use crate::BuildError; - -use super::compensated_sum; -use super::csr_matrix::CsrMatrix; -use super::factor::{factor_sparse, local_solver_build, ReducedFactor}; -use super::schur; - -/// Minimum number of rows to trigger parallel back-substitution. -const PAR_BACKSUB_THRESHOLD: usize = 10_000; -const PAR_BACKSUB_CHUNK: usize = 4096; - -#[inline] -fn subtract_mean(values: &mut [f64]) { - if values.is_empty() { - return; - } - let mean = compensated_sum(values) / values.len() as f64; - for val in values.iter_mut() { - *val -= mean; - } -} - -#[inline] -fn scale_by_diag_in_place(slice: &mut [f64], diagonal: &[f64]) { - debug_assert!(diagonal.len() >= slice.len()); - for (value, &scale) in slice.iter_mut().zip(diagonal.iter()) { - *value *= scale; - } -} - -/// Back-substitute for the eliminated block from a pre-scaled RHS. -fn backsub_block_from_scaled_rhs( - sol_output: &mut [f64], - scaled_rhs: &[f64], - cross_matrix: &CsrBlock, - inv_diag: &[f64], - sol_source: &[f64], - allow_inner_parallelism: bool, -) { - let n = sol_output.len(); - debug_assert!(scaled_rhs.len() >= n); - let recover = |i: usize| -> f64 { - let start = cross_matrix.indptr[i] as usize; - let end = cross_matrix.indptr[i + 1] as usize; - let mut sum = 0.0; - for idx in start..end { - let j = cross_matrix.indices[idx] as usize; - sum += cross_matrix.data[idx] * sol_source[j]; - } - scaled_rhs[i] + (inv_diag[i] * sum) - }; - if n > PAR_BACKSUB_THRESHOLD && allow_inner_parallelism { - sol_output - .par_chunks_mut(PAR_BACKSUB_CHUNK) - .enumerate() - .for_each(|(chunk_idx, chunk)| { - let row_start = chunk_idx * PAR_BACKSUB_CHUNK; - for (local_i, si) in chunk.iter_mut().enumerate() { - *si = recover(row_start + local_i); - } - }); - } else { - for (i, out) in sol_output.iter_mut().enumerate() { - *out = recover(i); - } - } -} - -/// Local subdomain solver using block elimination on the bipartite SDDM. -#[derive(Clone, serde::Serialize)] -pub struct BlockElimSolver { - /// Bipartite Gramian `C` and `Cᵀ`; diagonals are folded in at build time and not retained. - cross_tab: Arc, - /// `1 / D_elim[k]` for the eliminated (larger) diagonal block. - inv_diag_elim: Vec, - /// Reduced-system factor backend. - reduced_factor: ReducedFactor, - /// Internal DOF count; the matrix is single-sized, a cover living in `reduced_factor`. - #[serde(skip)] - n_internal: usize, - /// Internal factor dimension, including backend-added auxiliary vertices. - #[serde(skip)] - n_reduced: usize, - /// Original-to-SDDM coordinate map. - coordinates: CoordinateMap, -} - -impl<'de> serde::Deserialize<'de> for BlockElimSolver { - /// Validates every cross-field invariant [`Self::new`] assumes, against untrusted bytes. - fn deserialize>(deserializer: D) -> Result { - use serde::de::Error; - - #[derive(serde::Deserialize)] - struct Helper { - cross_tab: CrossTab, - inv_diag_elim: Vec, - reduced_factor: ReducedFactor, - coordinates: CoordinateMap, - } - - let h = Helper::deserialize(deserializer)?; - - // `ct`'s row count is the only witness bounding `c.ncols`. - let (c, ct) = h.cross_tab.into_parts(); - if !c.is_structurally_valid() { - return Err(D::Error::custom( - "cross_tab.c is not a structurally valid CSR block", - )); - } - if ct.nrows != c.ncols || ct.ncols != c.nrows || !ct.is_structurally_valid() { - return Err(D::Error::custom("cross_tab.ct shape disagrees with c")); - } - // Components are eliminated-major: the first block is scaled by `inv_diag_elim`. - let (n_eliminated, n_kept) = (c.nrows, c.ncols); - let n_internal = n_eliminated + n_kept; - - if h.inv_diag_elim.len() != n_eliminated { - return Err(D::Error::custom( - "inv_diag_elim length disagrees with eliminated block", - )); - } - if let CoordinateMap::Scaled(factors) = &h.coordinates { - if factors.len() != n_internal { - return Err(D::Error::custom( - "Scaled coordinate map length disagrees with n_internal", - )); - } - } - - if !h.reduced_factor.spans_kept_block(n_kept) { - return Err(D::Error::custom( - "reduced factor input dimension disagrees with the kept block", - )); - } - - // Rebuild from the validated `c` so the transpose cannot disagree. - Ok(BlockElimSolver::new( - CrossTab::eager(c), - h.inv_diag_elim, - h.reduced_factor, - h.coordinates, - )) - } -} - -/// An eliminated-major SDDM plus its eliminated-diagonal reciprocals, which it does not carry. -struct Eliminated { - matrix: SddmMatrix, - inv_diagonal: Vec, -} - -impl Eliminated { - fn new(matrix: SddmMatrix) -> Result { - debug_assert!( - matrix.n_eliminated() >= matrix.n_kept(), - "component is not eliminated-major" - ); - let inv_diagonal = matrix.diagonal[..matrix.n_eliminated()] - .iter() - .enumerate() - .map(|(i, &d)| { - if d > 0.0 { - Ok(1.0 / d) - } else { - Err(BuildError::SingularDiagonal { index: i }) - } - }) - .collect::, BuildError>>()?; - Ok(Self { - matrix, - inv_diagonal, - }) - } - - /// Fold of the signed matrix's Gremban cover; transient, dropped once its complement is built. - fn cover(&self) -> Result { - Self::new(assemble_bipartite_cover(&self.matrix)) - } - - /// A dense Cholesky spends nothing on sparsity, so it entails the exact complement. - fn factor_reduced( - fold: impl Borrow, - config: &LocalSolverConfig, - ) -> Result { - let this = fold.borrow(); - let exact_below = config.dense_threshold; - let exact = (exact_below > 0 && this.matrix.n_kept() <= exact_below) - .then(|| schur::exact_for_factor(&this.matrix, &this.inv_diagonal)); - if let Some(exact) = &exact { - match factor_complement(exact, config, ExactFailure::Error) { - Err(approx_chol::Error::DenseFactorizationFailed { .. }) => {} - result => return result.map_err(local_solver_build), - } - } - let complement = match &config.schur { - SchurMode::Approximate(cfg) => schur::sampled(&this.matrix, cfg), - SchurMode::Exact => { - exact.unwrap_or_else(|| schur::exact_for_factor(&this.matrix, &this.inv_diagonal)) - } - }; - // An owned fold is the transient cover; it is freed before the factor's fill is allocated. - drop(fold); - factor_complement(&complement, config, ExactFailure::FallBackToApproximate) - .map_err(local_solver_build) - } -} - -fn factor_complement( - complement: &CsrMatrix, - config: &LocalSolverConfig, - on_failure: ExactFailure, -) -> Result { - let approx_chol = config - .approx_chol - .to_approx_chol(config.dense_threshold, on_failure); - factor_sparse(complement, approx_chol) -} - -/// Gremban cover: SDDM, and acts on the antisymmetric `[z, -z]` subspace as the original. -fn assemble_bipartite_cover(matrix: &SddmMatrix) -> SddmMatrix { - let c = &matrix.cross_tab.c; - let n_rows = c.nrows; - let n_cols = c.ncols; - let n_cols_u32 = u32::try_from(n_cols).expect("cover columns exceed u32::MAX"); - // Both doubled sizes must fit u32: the cross-sheet shift emits columns up to `2*n_cols - 1`. - u32::try_from(2 * n_rows).expect("doubled cover rows exceed u32::MAX"); - u32::try_from(2 * n_cols).expect("doubled cover columns exceed u32::MAX"); - - let mut indptr = Vec::with_capacity(2 * n_rows + 1); - let mut indices = Vec::with_capacity(2 * c.nnz()); - let mut data = Vec::with_capacity(2 * c.nnz()); - indptr.push(0u32); - for copy_shifted in [false, true] { - for i in 0..n_rows { - let start = c.indptr[i] as usize; - let end = c.indptr[i + 1] as usize; - // Same-sheet columns precede cross-sheet ones so each output row stays column-sorted. - for column_shifted in [false, true] { - let column_base = if column_shifted { n_cols_u32 } else { 0 }; - let select_negative = column_shifted != copy_shifted; - for idx in start..end { - let value = c.data[idx]; - if (value < 0.0) != select_negative { - continue; - } - indices.push(c.indices[idx] + column_base); - data.push(value.abs()); - } - } - indptr.push(u32::try_from(indices.len()).expect("cover nonzeros exceed u32::MAX")); - } - } - let cover_c = CsrBlock { - indptr, - indices, - data, - nrows: 2 * n_rows, - ncols: 2 * n_cols, - }; - SddmMatrix { - cross_tab: CrossTab::new(cover_c), - diagonal: double_for_cover(&matrix.diagonal, n_rows), - ground_edges: double_for_cover(&matrix.ground_edges, n_rows), - grounding: matrix.grounding, - } -} - -/// Each Gremban sheet copies every vertex, so a per-vertex array doubles within each side. -fn double_for_cover(values: &[f64], n_rows: usize) -> Vec { - let (rows, cols) = values.split_at(n_rows); - [rows, rows, cols, cols].concat() -} - -impl BlockElimSolver { - /// Size of the block the elimination removes, the leading one. - fn n_eliminated(&self) -> usize { - self.cross_tab.n_rows() - } - - /// Size of the block the reduced factor solves. - fn n_kept(&self) -> usize { - self.cross_tab.n_cols() - } - - pub(crate) fn new( - cross_tab: impl Into>, - inv_diag_elim: Vec, - reduced_factor: ReducedFactor, - coordinates: CoordinateMap, - ) -> Self { - let cross_tab = cross_tab.into(); - // Every solve reads the transpose, so it is built here rather than inside the first. - cross_tab.ct(); - debug_assert!(reduced_factor.spans_kept_block(cross_tab.n_cols())); - debug_assert!( - reduced_factor.grounding() != Some(Grounding::Floating) - || reduced_factor.solve_dimension() == cross_tab.n_cols(), - "a floating factor has no ground slots for its mean to read", - ); - let n_internal = cross_tab.n_local(); - let n_reduced = reduced_factor.solve_dimension(); - Self { - cross_tab, - inv_diag_elim, - reduced_factor, - n_internal, - n_reduced, - coordinates, - } - } - - /// A signed component keeps only the cover's factor, keeping the matrix single-sized (#91). - pub(crate) fn build( - component: LocalComponent, - config: &LocalSolverConfig, - ) -> Result { - let LocalComponent { - matrix, - form, - coordinates, - } = component; - let eliminated = Eliminated::new(matrix)?; - - let factor = match form { - MatrixForm::Laplacian => { - let factor = ReducedFactor::Direct { - factor: Eliminated::factor_reduced(&eliminated, config)?, - grounding: eliminated.matrix.grounding, - }; - debug_assert!(factor.solve_dimension() >= factor.input_dimension()); - factor - } - // Surplus survives the cover, so it grounds as the signed matrix did. - MatrixForm::SignedPendingCover => ReducedFactor::Cover { - inner: Eliminated::factor_reduced(eliminated.cover()?, config)?, - m: eliminated.matrix.n_kept(), - }, - }; - - Ok(BlockElimSolver::new( - eliminated.matrix.cross_tab, - eliminated.inv_diagonal, - factor, - coordinates, - )) - } - - /// Components arrive oriented, so the row block is the eliminated side. - fn eliminate_and_recover( - &self, - rhs: &mut [f64], - sol: &mut [f64], - workspace: &mut [f64], - allow_inner_parallelism: bool, - ) -> Result<(), LocalSolveError> { - let (n_elim, n_keep) = (self.n_eliminated(), self.n_kept()); - let explicit_ground = self.reduced_factor.explicit_ground_index(n_keep); - debug_assert!( - n_keep <= self.n_reduced, - "reduced region must cover the kept block", - ); - debug_assert!( - explicit_ground.is_none_or(|g| n_keep <= g && g < self.n_reduced), - "grounded gauge slot must lie past the kept block inside the reduced system", - ); - - scale_by_diag_in_place(&mut rhs[..n_elim], &self.inv_diag_elim); - - let (reduced, embed) = workspace.split_at_mut(self.n_reduced); - self.cross_tab.ct().spmv_assign_add( - &rhs[..n_elim], - &rhs[n_elim..], - &mut reduced[..n_keep], - allow_inner_parallelism, - ); - match self.reduced_factor.grounding() { - Some(Grounding::Floating) => subtract_mean(reduced), - // The ground node's potential is the gauge subtracted after the solve. - Some(Grounding::Grounded) => { - if let Some(ground) = explicit_ground { - reduced[ground] = -compensated_sum(&reduced[..n_keep]); - } - } - // A cover grounds itself via the antisymmetric `[b, -b]` embed. - None => {} - } - - self.reduced_factor.solve_in_place(reduced, embed)?; - let (sol_elim, sol_kept) = sol.split_at_mut(n_elim); - sol_kept.copy_from_slice(&reduced[..n_keep]); - if let Some(ground) = explicit_ground { - let ground = reduced[ground]; - for v in sol_kept.iter_mut() { - *v -= ground; - } - } - - // Back-substitution reads the gauged kept block, so the gauge subtraction runs first. - backsub_block_from_scaled_rhs( - sol_elim, - &rhs[..n_elim], - &self.cross_tab.c, - &self.inv_diag_elim, - sol_kept, - allow_inner_parallelism, - ); - Ok(()) - } -} - -impl LocalSolver for BlockElimSolver { - fn n_local(&self) -> usize { - self.n_internal - } - - fn workspace_len(&self) -> usize { - self.n_internal + self.n_reduced + self.reduced_factor.embed_len() - } - - fn inner_parallelism_work_estimate(&self) -> usize { - let max_rows = self.cross_tab.n_rows().max(self.cross_tab.n_cols()); - if max_rows <= PAR_BACKSUB_THRESHOLD.max(PAR_SPMV_THRESHOLD) { - return 0; - } - - let cross_nnz = self.cross_tab.c.nnz(); - (2 * cross_nnz) + self.n_internal - } - - fn solve_local( - &self, - x: &mut [f64], - workspace: &mut [f64], - allow_inner_parallelism: bool, - ) -> Result<(), LocalSolveError> { - let n_eliminated = self.n_eliminated(); - // The elimination keeps its right-hand side apart from the solution it writes. - let (rhs, workspace) = workspace.split_at_mut(self.n_internal); - rhs.copy_from_slice(x); - let sol = x; - - self.coordinates.fold(rhs, n_eliminated); - if self.reduced_factor.grounding() == Some(Grounding::Floating) { - subtract_mean(rhs); - } - - self.eliminate_and_recover(rhs, sol, workspace, allow_inner_parallelism)?; - - if self.reduced_factor.grounding() == Some(Grounding::Floating) { - subtract_mean(sol); - } - self.coordinates.unfold(sol, n_eliminated); - Ok(()) - } -} - -#[cfg(test)] -mod tests; diff --git a/crates/within/src/block_elim/solver/tests.rs b/crates/within/src/block_elim/solver/tests.rs deleted file mode 100644 index 1a4bcc74..00000000 --- a/crates/within/src/block_elim/solver/tests.rs +++ /dev/null @@ -1,442 +0,0 @@ -use rstest::rstest; - -use super::*; - -use crate::block_elim::csr_matrix::CsrMatrix; -use crate::config::{ApproxCholConfig, SchurMode, DEFAULT_DENSE_SCHUR_THRESHOLD}; - -impl BlockElimSolver { - /// Solves `r` with the workspace preset to `stale`. - pub(crate) fn solve_for_test(&self, r: &[f64], stale: f64) -> Vec { - let mut x = r.to_vec(); - let mut workspace = vec![stale; self.workspace_len()]; - self.solve_local(&mut x, &mut workspace, false) - .expect("solve_local failed"); - x - } -} - -#[test] -fn test_subtract_mean_empty() { - let mut data: Vec = Vec::new(); - subtract_mean(&mut data); - assert!(data.is_empty()); -} - -#[test] -fn test_subtract_mean_basic() { - let mut data = vec![2.0, 4.0, 6.0]; - subtract_mean(&mut data); - // mean = 4.0 - assert!((data[0] - (-2.0)).abs() < 1e-14); - assert!((data[1] - 0.0).abs() < 1e-14); - assert!((data[2] - 2.0).abs() < 1e-14); -} - -#[test] -fn fold_inverts_the_eliminated_diagonal() { - // No cross entries, so the fold reads the diagonal alone. - let matrix = SddmMatrix::from_dense_for_test( - &[0.0; 6], - 3, - 2, - vec![5.0, 6.0, 8.0, 7.0, 9.0], - Grounding::Grounded, - ); - let eliminated = Eliminated::new(matrix).expect("a positive diagonal must fold"); - - assert_eq!(eliminated.inv_diagonal.len(), 3); - for (&got, &expected) in eliminated - .inv_diagonal - .iter() - .zip([1.0 / 5.0, 1.0 / 6.0, 1.0 / 8.0].iter()) - { - assert!((got - expected).abs() < 1e-12); - } -} - -#[test] -fn fold_rejects_a_zero_eliminated_diagonal() { - let matrix = SddmMatrix::from_dense_for_test( - &[0.0; 6], - 3, - 2, - vec![5.0, 6.0, 0.0, 8.0, 9.0], - Grounding::Grounded, - ); - - match Eliminated::new(matrix) { - Err(BuildError::SingularDiagonal { index: 2, .. }) => {} - Err(e) => panic!("expected SingularDiagonal at index 2, got: {e}"), - Ok(_) => panic!("expected SingularDiagonal error, got Ok"), - } -} - -#[rstest] -fn an_unusable_dense_pivot_is_retried_rather_than_fatal( - #[values(SchurMode::Exact, SchurMode::Approximate(Default::default()))] schur: SchurMode, -) { - // The exact attempt only fails when weight spread costs the complement its definiteness. - let matrix = - SddmMatrix::laplacian_for_test(&[0.0, 1e-16, 1e2, 1e0, 1e5, 0.0, 0.0, 1e14, 0.0], 3, 3); - let eliminated = Eliminated::new(matrix).expect("a positive diagonal must fold"); - let exact = schur::exact_for_factor(&eliminated.matrix, &eliminated.inv_diagonal); - let exact_only = ApproxCholConfig::default() - .to_approx_chol(DEFAULT_DENSE_SCHUR_THRESHOLD, ExactFailure::Error); - assert!( - matches!( - factor_sparse(&exact, exact_only), - Err(approx_chol::Error::DenseFactorizationFailed { .. }) - ), - "the fixture no longer reaches the fall-through" - ); - - let config = LocalSolverConfig { - approx_chol: ApproxCholConfig::default(), - schur, - dense_threshold: DEFAULT_DENSE_SCHUR_THRESHOLD, - ..Default::default() - }; - assert!( - Eliminated::factor_reduced(&eliminated, &config).is_ok(), - "an unusable pivot must not be fatal" - ); -} - -/// Build an eliminated-major CrossTab (`n_rows > n_cols`, as orientation -/// guarantees) whose two cross entries leave three isolated rows, plus its -/// build-time diagonal. -fn make_cross_tab() -> (CrossTab, Vec) { - let c_dense = vec![ - 1.0, 0.0, // row 0 - 0.0, 1.0, // row 1 - 0.0, 0.0, // row 2 - 0.0, 0.0, // row 3 - 0.0, 0.0, // row 4 - ]; - let diagonal = vec![2.0, 3.0, 1.0, 1.0, 1.0, 2.0, 3.0]; - (CrossTab::from_dense_for_test(&c_dense, 5, 2), diagonal) -} - -#[test] -fn block_elim_solver_solves_a_two_block_component() { - let (cross_tab, diagonals) = make_cross_tab(); - assert_eq!(cross_tab.n_rows(), 5); - assert_eq!(cross_tab.n_cols(), 2); - - let config = LocalSolverConfig { - approx_chol: ApproxCholConfig::default(), - schur: SchurMode::Exact, - dense_threshold: 0, // disable dense fast path to ensure sparse path is covered - ..Default::default() - }; - let component = LocalComponent::general_for_test(cross_tab, diagonals); - let solver = BlockElimSolver::build(component, &config).expect("block-elim build failed"); - - let n_local = solver.n_local(); - assert_eq!(n_local, 7); - - let r: Vec = (0..n_local).map(|i| (i as f64 + 1.0) * 0.5).collect(); - let sol = solver.solve_for_test(&r, 0.0); - - for (i, &v) in sol.iter().enumerate() { - assert!(v.is_finite(), "sol[{i}] = {v} is not finite"); - } - let sol_norm: f64 = sol.iter().map(|v| v * v).sum::().sqrt(); - assert!(sol_norm > 1e-15, "solution is unexpectedly all-zero"); -} - -#[test] -fn grounded_two_block_solve_is_leak_free() { - // Grounded sparse factor: pins no-leak (pre-dirtying must not move the result). - let (cross_tab, diagonals) = make_cross_tab(); - let config = LocalSolverConfig { - approx_chol: ApproxCholConfig::default(), - schur: SchurMode::Exact, - dense_threshold: 0, // force the sparse grounded path → explicit ground vertex - ..Default::default() - }; - let component = LocalComponent::general_for_test(cross_tab, diagonals); - let solver = BlockElimSolver::build(component, &config).expect("block-elim build failed"); - // Original bipartite Gram A the solver inverts (diagonals + the two cross entries). - let n = solver.n_local(); - let mut a = vec![vec![0.0; n]; n]; - for (i, d) in [2.0, 3.0, 1.0, 1.0, 1.0, 2.0, 3.0].into_iter().enumerate() { - a[i][i] = d; - } - for (i, j) in [(0, 5), (1, 6)] { - a[i][j] = 1.0; - a[j][i] = 1.0; - } - - let r = [0.5, 3.0, -1.25, 0.75, 2.0, 1.0, -2.0]; - let clean = solver.solve_for_test(&r, 0.0); - let dirty = solver.solve_for_test(&r, 1e6); - assert_eq!( - clean, dirty, - "pre-dirtied overlap/scratch slots leaked into the solution", - ); - - for i in 0..n { - let ax: f64 = (0..n).map(|j| a[i][j] * clean[j]).sum(); - assert!( - (ax - r[i]).abs() < 1e-9, - "row {i}: A·x = {ax}, expected {}", - r[i], - ); - } -} - -/// Live 1×1 components keep n_keep = 0, so the solve degenerates to x = r/d exactly. -#[rstest] -#[case::row_singleton(1, 0)] -#[case::column_singleton(0, 1)] -fn trivial_singleton_component_solves_r_over_d(#[case] n_rows: usize, #[case] n_cols: usize) { - let config = LocalSolverConfig { - approx_chol: ApproxCholConfig::default(), - schur: SchurMode::Exact, - dense_threshold: 0, - ..Default::default() - }; - let cross_tab = CrossTab::from_dense_for_test(&[], n_rows, n_cols); - let diagonal = vec![4.0; n_rows + n_cols]; - let component = LocalComponent::general_for_test(cross_tab, diagonal); - let solver = BlockElimSolver::build(component, &config).expect("trivial 1×1 build"); - assert_eq!(solver.n_local(), 1); - - let sol = solver.solve_for_test(&[2.0], 0.0); - assert_eq!(sol[0], 0.5, "expected r/d"); -} - -#[test] -fn sampled_sparse_preserves_barely_pd_direction() { - let surplus = 5e-10; - let cross_tab = CrossTab::from_dense_for_test(&[1.0], 1, 1); - let component = LocalComponent::general_for_test(cross_tab, vec![1.0 + surplus, 1.0]); - let config = LocalSolverConfig { - approx_chol: ApproxCholConfig::default(), - schur: SchurMode::Approximate(crate::config::ApproxSchurConfig::default()), - dense_threshold: 0, - ..Default::default() - }; - let solver = BlockElimSolver::build(component, &config).unwrap(); - let solution = solver.solve_for_test(&[1.0, -1.0], 0.0); - - let ax = [ - (1.0 + surplus) * solution[0] + solution[1], - solution[0] + solution[1], - ]; - assert!((ax[0] - 1.0).abs() < 1e-6); - assert!((ax[1] + 1.0).abs() < 1e-6); -} - -fn two_dof_sddm(diagonal: [f64; 2], off_diagonal: f64) -> CsrMatrix { - let data = vec![diagonal[0], off_diagonal, off_diagonal, diagonal[1]]; - CsrMatrix::new(vec![0, 2, 4], vec![0, 1, 0, 1], data, 2) -} - -fn explicit_ground_laplacian() -> CsrMatrix { - let (large, small) = (1e7, 1e-9); - let principal = two_dof_sddm([large, large + small], -large); - schur::build_explicit_laplacian(&principal, &[0.0, small], Grounding::Grounded) -} - -#[rstest] -// approx-chol declines to ground a surplus below its resolvable pivot scale. -#[case::explicit_ground(explicit_ground_laplacian(), 3)] -// A resolvable surplus makes approx-chol append a ground vertex of its own. -#[case::backend_ground(two_dof_sddm([2.5, 2.5], -2.0), 2)] -fn stale_ground_slots_past_the_kept_block_never_reach_the_solution( - #[case] reduced: CsrMatrix, - #[case] input_dimension: usize, - #[values(3.0, -7.0)] stale: f64, -) { - let config = ApproxCholConfig::default().to_approx_chol( - DEFAULT_DENSE_SCHUR_THRESHOLD, - ExactFailure::FallBackToApproximate, - ); - let factor = ReducedFactor::Direct { - factor: factor_sparse(&reduced, config).expect("factorization must succeed"), - grounding: Grounding::Grounded, - }; - assert_eq!( - (factor.input_dimension(), factor.solve_dimension()), - (input_dimension, 3) - ); - - let cross_tab = CrossTab::from_dense_for_test(&[0.0; 6], 3, 2); - let solver = BlockElimSolver::new(cross_tab, vec![1.0; 3], factor, CoordinateMap::default()); - - let r = [0.0, 0.0, 0.0, 1.0, -1.0]; - assert_eq!( - solver.solve_for_test(&r, stale), - solver.solve_for_test(&r, 0.0) - ); -} - -/// Stars of ≤2 entries sample deterministically-exact, so all three arms share one oracle. -#[rstest] -#[case::dense_full_minor(8, SchurMode::Exact)] -#[case::exact_sparse(0, SchurMode::Exact)] -#[case::sampled_sparse(0, SchurMode::Approximate(crate::config::ApproxSchurConfig::default()))] -fn signed_component_realizes_congruence_transformed_solve( - #[case] dense_threshold: usize, - #[case] schur: SchurMode, -) { - let (n_rows, n_cols) = (3usize, 2usize); - let d: Vec = vec![-1.0, 4.0, 0.25, 2.0, -0.5]; - let c_hat = [[1.0, 3.0], [2.0, 0.0], [0.5, 1.5]]; - let diag_hat = [4.0, 2.5, 2.0, 4.0, 5.0]; // surplus on q0, q2, r0 - - let mut c_raw = vec![0.0; n_rows * n_cols]; - for i in 0..n_rows { - for j in 0..n_cols { - c_raw[i * n_cols + j] = c_hat[i][j] / (d[i] * d[n_rows + j]); - } - } - - let cross_tab = CrossTab::from_dense_for_test(&c_raw, n_rows, n_cols); - let diagonals: Vec = (0..n_rows + n_cols) - .map(|k| diag_hat[k] / (d[k] * d[k])) - .collect(); - let config = LocalSolverConfig { - approx_chol: ApproxCholConfig::default(), - schur, - dense_threshold, - ..Default::default() - }; - // SDDM factors fold the bipartite negation in: f = d on q, -d on r. - let factors: Vec = d - .iter() - .enumerate() - .map(|(i, &v)| if i < n_rows { v } else { -v }) - .collect(); - let component = LocalComponent::with_factors_for_test(cross_tab, diagonals, &factors); - let solver = - BlockElimSolver::build(component, &config).expect("signed block-elim build failed"); - - let n = n_rows + n_cols; - let r = [0.5, 3.0, -1.25, 1.0, -2.0]; - let sol = solver.solve_for_test(&r, 0.0); - - for i in 0..n { - let mut ax = 0.0; - for j in 0..n { - let a_ij = if i == j { - diag_hat[i] / (d[i] * d[i]) - } else if i < n_rows && j >= n_rows { - c_raw[i * n_cols + (j - n_rows)] - } else if j < n_rows && i >= n_rows { - c_raw[j * n_cols + (i - n_rows)] - } else { - 0.0 - }; - ax += a_ij * sol[j]; - } - assert!( - (ax - r[i]).abs() < 1e-9, - "row {i}: A·x = {ax}, expected {}", - r[i] - ); - } -} - -#[test] -fn frustrated_component_solves_exactly_through_cover() { - // Dense only: it is the one arm whose factor is exact, so the oracle A·x = r is sharp. - let (n_rows, n_cols) = (2usize, 2usize); - let c_raw = [1.0, 1.5, 2.0, -1.0]; - let a = [ - [2.5, 0.0, 1.0, 1.5], - [0.0, 3.0, 2.0, -1.0], - [1.0, 2.0, 3.4, 0.0], - [1.5, -1.0, 0.0, 2.9], - ]; - - let component = LocalComponent::general_for_test( - CrossTab::from_dense_for_test(&c_raw, n_rows, n_cols), - vec![a[0][0], a[1][1], a[2][2], a[3][3]], - ); - let config = LocalSolverConfig { - approx_chol: ApproxCholConfig::default(), - schur: SchurMode::Exact, - dense_threshold: 8, - ..Default::default() - }; - let solver = - BlockElimSolver::build(component, &config).expect("covered block-elim build failed"); - let n = n_rows + n_cols; - // #91: the stored operator stays single-sized; the cover lives inside the reduced factor. - assert_eq!(solver.n_local(), n, "operator stays single-sized"); - assert_eq!( - solver.cross_tab.n_local(), - n, - "stored cross-tab not doubled" - ); - assert!( - matches!(solver.reduced_factor, ReducedFactor::Cover { .. }), - "frustrated component must reduce through a cover", - ); - - let r = [1.0, -2.0, 0.5, 3.0]; - let sol = solver.solve_for_test(&r, 0.0); - - for i in 0..n { - let ax: f64 = (0..n).map(|j| a[i][j] * sol[j]).sum(); - assert!( - (ax - r[i]).abs() < 1e-9, - "row {i}: A·x = {ax}, expected {}", - r[i] - ); - } -} - -/// A structurally valid solver built the normal way; the deserialize-validation -/// tests below tamper one field of a clone and assert the round-trip rejects it. -fn valid_solver_for_deser() -> BlockElimSolver { - let (cross_tab, diagonals) = make_cross_tab(); - let config = LocalSolverConfig { - approx_chol: ApproxCholConfig::default(), - schur: SchurMode::Exact, - dense_threshold: 0, - ..Default::default() - }; - let component = LocalComponent::general_for_test(cross_tab, diagonals); - BlockElimSolver::build(component, &config).expect("block-elim build failed") -} - -#[test] -fn valid_solver_round_trips() { - let solver = valid_solver_for_deser(); - let bytes = postcard::to_stdvec(&solver).expect("serialize"); - let restored: BlockElimSolver = postcard::from_bytes(&bytes).expect("deserialize"); - assert_eq!(restored.n_local(), solver.n_local()); - assert_eq!(restored.workspace_len(), solver.workspace_len()); -} - -#[test] -fn inv_diag_elim_length_mismatch_is_rejected() { - let mut bad = valid_solver_for_deser(); - bad.inv_diag_elim.push(0.0); - let bytes = postcard::to_stdvec(&bad).expect("serialize"); - assert!(postcard::from_bytes::(&bytes).is_err()); -} - -#[test] -fn scaled_coordinate_length_mismatch_is_rejected() { - let mut bad = valid_solver_for_deser(); - bad.coordinates = CoordinateMap::Scaled(vec![1.0; bad.n_internal + 3].into_boxed_slice()); - let bytes = postcard::to_stdvec(&bad).expect("serialize"); - assert!(postcard::from_bytes::(&bytes).is_err()); -} - -#[test] -fn kept_block_wider_than_the_reduced_factor_is_rejected() { - let mut bad = valid_solver_for_deser(); - // Earlier cross-field witnesses stay intact, so the factor's span is the check under test. - let mut c = bad.cross_tab.c.clone(); - c.ncols += 2; - bad.cross_tab = Arc::new(CrossTab::eager(c)); - let bytes = postcard::to_stdvec(&bad).expect("serialize"); - assert!(postcard::from_bytes::(&bytes).is_err()); -} diff --git a/crates/within/src/block_elim/tests.rs b/crates/within/src/block_elim/tests.rs deleted file mode 100644 index fa9c8d30..00000000 --- a/crates/within/src/block_elim/tests.rs +++ /dev/null @@ -1,34 +0,0 @@ -use super::compensated_sum; - -// A naive left-to-right fold drops the small middle term (`1e16 + 1.0` rounds -// back to `1e16`), returning 0.0; Neumaier compensation recovers the dropped -// bits and returns exactly 1.0. The two orderings exercise both dominance -// branches of the running compensation (`sum.abs() >= value.abs()` and `else`), -// so each carried term is pinned rather than the naive and compensated results -// coinciding. -#[test] -fn recovers_catastrophic_cancellation_in_both_branches() { - // Large value first: the `sum.abs() >= value.abs()` branch carries the 1.0. - let big_first = [1e16, 1.0, -1e16]; - assert_eq!( - big_first.iter().sum::(), - 0.0, - "naive fold loses the term" - ); - assert_eq!(compensated_sum(&big_first), 1.0); - - // Small value first: the `else` branch carries the 1.0. - let small_first = [1.0, 1e16, -1e16]; - assert_eq!( - small_first.iter().sum::(), - 0.0, - "naive fold loses the term" - ); - assert_eq!(compensated_sum(&small_first), 1.0); -} - -#[test] -fn handles_empty_and_singleton() { - assert_eq!(compensated_sum(&[]), 0.0); - assert_eq!(compensated_sum(&[42.0]), 42.0); -} diff --git a/crates/within/src/block_solver.rs b/crates/within/src/block_solver.rs new file mode 100644 index 00000000..8b06fe80 --- /dev/null +++ b/crates/within/src/block_solver.rs @@ -0,0 +1,60 @@ +//! Schwarz local solver over approx-chol's irreducible-block factors. + +use approx_chol::{BlockFactor, InnerParallelism}; +use schwarz_precond::{LocalSolveError, LocalSolver}; +use serde::{Deserialize, Serialize}; + +/// Consecutive irreducible blocks; the subdomain's local order concatenates their vertices. +#[derive(Clone, Serialize, Deserialize)] +pub(crate) struct BlockSolver { + blocks: Vec, +} + +impl BlockSolver { + pub(crate) fn new(blocks: Vec) -> Self { + Self { blocks } + } +} + +impl LocalSolver for BlockSolver { + fn n_local(&self) -> usize { + self.blocks.iter().map(|b| b.vertices().len()).sum() + } + + fn workspace_len(&self) -> usize { + self.blocks + .iter() + .map(BlockFactor::scratch_len) + .max() + .unwrap_or(0) + } + + fn solve_local( + &self, + x: &mut [f64], + workspace: &mut [f64], + allow_inner_parallelism: bool, + ) -> Result<(), LocalSolveError> { + let inner = if allow_inner_parallelism { + InnerParallelism::Allow + } else { + InnerParallelism::Forbid + }; + let mut rest = x; + for block in &self.blocks { + let (local, tail) = rest.split_at_mut(block.vertices().len()); + block + .solve_local(local, &mut workspace[..block.scratch_len()], inner) + .map_err(|e| LocalSolveError::BackendFailed { + context: "BlockSolver::solve_local", + message: e.to_string(), + })?; + rest = tail; + } + Ok(()) + } + + fn inner_parallelism_work_estimate(&self) -> usize { + self.blocks.iter().map(BlockFactor::splittable_work).sum() + } +} diff --git a/crates/within/src/config.rs b/crates/within/src/config.rs index 6a1db8b1..a2684b79 100644 --- a/crates/within/src/config.rs +++ b/crates/within/src/config.rs @@ -13,7 +13,7 @@ use crate::error::BuildError; pub use schwarz_precond::{Staleness, StalenessError}; -/// Default `n_keep` threshold below which a Schur domain tries the exact dense backend. +/// Default block dimension at or below which a block is factored exactly. pub(crate) const DEFAULT_DENSE_SCHUR_THRESHOLD: usize = 24; /// Configuration for approximate Cholesky factorization. @@ -25,24 +25,6 @@ pub struct ApproxCholConfig { pub split_merge: Option, } -impl ApproxCholConfig { - pub(crate) fn to_approx_chol( - self, - exact_below: usize, - on_failure: approx_chol::ExactFailure, - ) -> approx_chol::Config { - approx_chol::Config { - seed: self.seed, - split_merge: self.split_merge, - // `max_dim: 0` claims no block, which is what `dense_threshold: 0` asks for. - backend: approx_chol::Backend::ExactBelow { - max_dim: exact_below, - on_failure, - }, - } - } -} - /// Schur reduction mode: sampled bounds per-subdomain cost, exact trades speed for fidelity. #[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] pub enum SchurMode { @@ -65,7 +47,7 @@ pub struct LocalSolverConfig { pub approx_chol: ApproxCholConfig, /// Schur-complement reduction mode (default: approximate). pub schur: SchurMode, - /// Reduced-size threshold at or below which fill-in is affordable; `0` disables. + /// Largest Cholesky dimension factored exactly, a signed block's two cover sheets counted; `0` disables. pub dense_threshold: usize, /// Certification policy for the diagonal scaling of signed components. pub scaling: ScalingConfig, @@ -74,11 +56,33 @@ pub struct LocalSolverConfig { impl LocalSolverConfig { /// Separate from any build so a deferred escalation can reject the config at construction. pub(crate) fn validate(&self) -> Result<(), BuildError> { - let tolerance = self.scaling.tolerance; - if !tolerance.is_finite() || tolerance < 0.0 { - return Err(BuildError::InvalidScalingTolerance { value: tolerance }); + self.scaling_policy().map(|_| ()) + } + + pub(crate) fn scaling_policy(&self) -> Result { + let ScalingConfig { + tolerance, + max_iterations, + .. + } = self.scaling; + approx_chol::ScalingPolicy::new(tolerance, max_iterations) + .ok_or(BuildError::InvalidScalingTolerance { value: tolerance }) + } + + pub(crate) fn factor_config(&self) -> approx_chol::Config { + approx_chol::Config { + seed: self.approx_chol.seed, + split_merge: self.approx_chol.split_merge, + exact_max_dim: self.dense_threshold, + complement: match self.schur { + SchurMode::Approximate(ApproxSchurConfig { split }) => { + approx_chol::SchurMode::Approximate { + split_merge: Some(split), + } + } + SchurMode::Exact => approx_chol::SchurMode::Exact, + }, } - Ok(()) } } @@ -99,9 +103,9 @@ impl Default for LocalSolverConfig { /// Frustration is always a hard error; this governs only the dominance certificate. #[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] pub struct ScalingConfig { - /// Relative slack for weak dominance; PSD-boundary designs hover at ≈ 1e-12. + /// Relative slack for weak dominance, in `[0, 1)`; PSD-boundary designs hover at ≈ 1e-12. pub tolerance: f64, - /// Reduced-CG iteration budget; each iteration costs 2–4 matrix-vector passes. + /// Scaling iteration budget per irreducible block. pub max_iterations: usize, /// Disposition when certification fails. pub on_failure: ScalingFailure, @@ -126,18 +130,16 @@ impl Default for ScalingConfig { } } -/// Approximate Schur via GKS 2023 Algorithm 3: unbiased edge weights at O(deg) fill. +/// Approximate Schur via GKS 2023 Algorithm 3, drawing from [`ApproxCholConfig::seed`]. #[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)] pub struct ApproxSchurConfig { - /// Random seed for the clique-tree sampler. - pub seed: u64, /// Star-edge split factor: `1` is standard, `k > 1` a denser approximation at more fill-in. pub split: u32, } impl Default for ApproxSchurConfig { fn default() -> Self { - Self { seed: 0, split: 1 } + Self { split: 1 } } } diff --git a/crates/within/src/csr_block.rs b/crates/within/src/csr_block.rs deleted file mode 100644 index 74efb178..00000000 --- a/crates/within/src/csr_block.rs +++ /dev/null @@ -1,368 +0,0 @@ -use rayon::prelude::*; - -/// Minimum number of rows to trigger parallel SpMV. -pub(crate) const PAR_SPMV_THRESHOLD: usize = 10_000; -/// Target number of non-zeros per parallel chunk. -const TARGET_NNZ_PER_CHUNK: usize = 32_768; - -/// Checked `usize -> u32`: a silent `as` truncation would corrupt the structure undiagnosed. -#[inline] -pub(crate) fn to_u32(x: usize) -> u32 { - u32::try_from(x).expect("CSR index exceeds u32::MAX") -} - -/// Rectangular CSR block storing `C` or `Cᵀ`, column indices ascending within each row. -#[derive(Clone, serde::Serialize, serde::Deserialize)] -pub(crate) struct CsrBlock { - pub(crate) indptr: Vec, - pub(crate) indices: Vec, - pub(crate) data: Vec, - pub(crate) nrows: usize, - pub(crate) ncols: usize, -} - -impl CsrBlock { - /// Number of stored non-zeros. - pub(crate) fn nnz(&self) -> usize { - self.data.len() - } - - /// `nrows`/`ncols` sit beside the arrays, so untrusted bytes can make them disagree. - pub(crate) fn is_structurally_valid(&self) -> bool { - if self.nrows.checked_add(1) != Some(self.indptr.len()) { - return false; - } - let nnz = self.indices.len(); - if self.data.len() != nnz { - return false; - } - if self.indptr[0] != 0 || *self.indptr.last().unwrap() as usize != nnz { - return false; - } - if self.indptr.windows(2).any(|w| w[0] > w[1]) { - return false; - } - self.indices.iter().all(|&j| (j as usize) < self.ncols) - } - - /// Row `i` as `(column, weight)` pairs. - pub(crate) fn row(&self, i: usize) -> impl Iterator + '_ { - let start = self.indptr[i] as usize; - let end = self.indptr[i + 1] as usize; - self.indices[start..end] - .iter() - .zip(&self.data[start..end]) - .map(|(&j, &w)| (j as usize, w)) - } - - /// Transpose in O(nnz); output rows come out sorted because source rows go ascending. - pub(crate) fn transpose(&self) -> CsrBlock { - let nnz = self.nnz(); - let mut row_counts = vec![0u32; self.ncols]; - for &col in &self.indices { - row_counts[col as usize] += 1; - } - let mut indptr = vec![0u32; self.ncols + 1]; - for i in 0..self.ncols { - indptr[i + 1] = indptr[i] + row_counts[i]; - } - let mut cursor = indptr[..self.ncols].to_vec(); - let mut indices = vec![0u32; nnz]; - let mut data = vec![0.0f64; nnz]; - for src_row in 0..self.nrows { - let start = self.indptr[src_row] as usize; - let end = self.indptr[src_row + 1] as usize; - for idx in start..end { - let dst_row = self.indices[idx] as usize; - let pos = cursor[dst_row] as usize; - indices[pos] = src_row as u32; - data[pos] = self.data[idx]; - cursor[dst_row] += 1; - } - } - CsrBlock { - indptr, - indices, - data, - nrows: self.ncols, - ncols: self.nrows, - } - } - - /// Build a CSR block from a row-major dense table (`table[i * ncols + j]`), skipping zeros. - pub(crate) fn from_dense_table(table: &[f64], nrows: usize, ncols: usize) -> Self { - debug_assert_eq!(table.len(), nrows * ncols); - let mut indptr = vec![0u32; nrows + 1]; - for i in 0..nrows { - let row_start = i * ncols; - let mut count = 0u32; - for j in 0..ncols { - if table[row_start + j] != 0.0 { - count += 1; - } - } - indptr[i + 1] = indptr[i] + count; - } - let nnz = indptr[nrows] as usize; - let mut indices = Vec::with_capacity(nnz); - let mut data = Vec::with_capacity(nnz); - for i in 0..nrows { - let row_start = i * ncols; - for j in 0..ncols { - let v = table[row_start + j]; - if v != 0.0 { - indices.push(j as u32); - data.push(v); - } - } - } - CsrBlock { - indptr, - indices, - data, - nrows, - ncols, - } - } - - /// `y = base + A x`, fusing the base copy with accumulation to avoid an extra pass. - pub(crate) fn spmv_assign_add( - &self, - x: &[f64], - base: &[f64], - y: &mut [f64], - allow_inner_parallelism: bool, - ) { - debug_assert!(x.len() >= self.ncols); - debug_assert!(base.len() >= self.nrows); - debug_assert!(y.len() >= self.nrows); - if self.nrows > PAR_SPMV_THRESHOLD && allow_inner_parallelism { - self.par_spmv_assign_add(x, base, y); - } else { - self.seq_spmv_assign_add(x, base, y); - } - } - - fn seq_spmv_assign_add(&self, x: &[f64], base: &[f64], y: &mut [f64]) { - for (i, yi) in y[..self.nrows].iter_mut().enumerate() { - let start = self.indptr[i] as usize; - let end = self.indptr[i + 1] as usize; - let row_data = &self.data[start..end]; - let row_idx = &self.indices[start..end]; - let mut acc = base[i]; - for (&val, &col) in row_data.iter().zip(row_idx) { - acc += val * x[col as usize]; - } - *yi = acc; - } - } - - fn par_spmv_assign_add(&self, x: &[f64], base: &[f64], y: &mut [f64]) { - let indptr = &self.indptr; - let indices = &self.indices; - let data = &self.data; - let nnz = self.nnz(); - let nrows = self.nrows; - let avg_nnz_per_row = nnz / nrows.max(1); - let chunk = (TARGET_NNZ_PER_CHUNK / avg_nnz_per_row.max(1)).clamp(256, 8192); - - y[..self.nrows] - .par_chunks_mut(chunk) - .enumerate() - .for_each(|(chunk_idx, y_chunk)| { - let row_start = chunk_idx * chunk; - for (local_i, yi) in y_chunk.iter_mut().enumerate() { - let i = row_start + local_i; - let start = indptr[i] as usize; - let end = indptr[i + 1] as usize; - let row_data = &data[start..end]; - let row_idx = &indices[start..end]; - let mut acc = base[i]; - for (&val, &col) in row_data.iter().zip(row_idx) { - acc += val * x[col as usize]; - } - *yi = acc; - } - }); - } -} - -#[cfg(test)] -mod tests { - use super::*; - - fn make_3x4_block() -> CsrBlock { - CsrBlock::from_dense_table( - &[1.0, 0.0, 2.0, 0.0, 0.0, 3.0, 0.0, 4.0, 5.0, 0.0, 6.0, 0.0], - 3, - 4, - ) - } - - #[test] - fn test_from_dense_table_structure() { - let b = make_3x4_block(); - assert_eq!(b.nrows, 3); - assert_eq!(b.ncols, 4); - assert_eq!(b.nnz(), 6); - - // Row 0: cols 0, 2 - assert_eq!(b.indptr[0], 0); - assert_eq!(b.indptr[1], 2); - assert_eq!(b.indices[0], 0); - assert_eq!(b.indices[1], 2); - assert_eq!(b.data[0], 1.0); - assert_eq!(b.data[1], 2.0); - - // Row 1: cols 1, 3 - assert_eq!(b.indptr[2], 4); - assert_eq!(b.indices[2], 1); - assert_eq!(b.indices[3], 3); - assert_eq!(b.data[2], 3.0); - assert_eq!(b.data[3], 4.0); - - // Row 2: cols 0, 2 - assert_eq!(b.indptr[3], 6); - assert_eq!(b.indices[4], 0); - assert_eq!(b.indices[5], 2); - assert_eq!(b.data[4], 5.0); - assert_eq!(b.data[5], 6.0); - } - - #[test] - fn test_from_dense_table_all_zeros() { - let b = CsrBlock::from_dense_table(&[0.0; 6], 2, 3); - assert_eq!(b.nrows, 2); - assert_eq!(b.ncols, 3); - assert_eq!(b.nnz(), 0); - assert_eq!(b.indptr, vec![0, 0, 0]); - } - - #[test] - fn test_transpose_basic() { - let b = make_3x4_block(); - let bt = b.transpose(); - - assert_eq!(bt.nrows, 4); - assert_eq!(bt.ncols, 3); - assert_eq!(bt.nnz(), 6); - - // A^T[j, i] == A[i, j] on the values, not just the structure. - - // Row 0 of transpose: columns 0 and 2 (from original rows 0 and 2 having col 0) - let r0_start = bt.indptr[0] as usize; - let r0_end = bt.indptr[1] as usize; - let r0_cols: Vec = bt.indices[r0_start..r0_end].to_vec(); - let r0_vals: Vec = bt.data[r0_start..r0_end].to_vec(); - assert_eq!(r0_cols, vec![0, 2]); // rows 0, 2 of original - assert_eq!(r0_vals, vec![1.0, 5.0]); // values at (0,0) and (2,0) of original - } - - #[test] - fn test_transpose_transpose_roundtrip() { - let b = make_3x4_block(); - let btt = b.transpose().transpose(); - - assert_eq!(btt.nrows, b.nrows); - assert_eq!(btt.ncols, b.ncols); - assert_eq!(btt.nnz(), b.nnz()); - assert_eq!(btt.indptr, b.indptr); - assert_eq!(btt.indices, b.indices); - // Values should match - for (a, e) in btt.data.iter().zip(b.data.iter()) { - assert!((a - e).abs() < 1e-14); - } - } - - #[test] - fn test_nnz() { - let b = make_3x4_block(); - assert_eq!(b.nnz(), 6); - - let empty = CsrBlock::from_dense_table(&[0.0; 4], 2, 2); - assert_eq!(empty.nnz(), 0); - } - - #[test] - fn test_spmv_assign_add_basic() { - let b = make_3x4_block(); - let x = vec![1.0, 1.0, 1.0, 1.0]; - let base = vec![10.0, 20.0, 30.0]; - let mut y = vec![0.0; 3]; - - b.spmv_assign_add(&x, &base, &mut y, true); - - assert!((y[0] - 13.0).abs() < 1e-14); - assert!((y[1] - 27.0).abs() < 1e-14); - assert!((y[2] - 41.0).abs() < 1e-14); - } - - #[test] - fn test_spmv_assign_add_zero_x_preserves_base() { - let b = make_3x4_block(); - let x = vec![0.0; 4]; - let base = vec![5.0; 3]; - let mut y = vec![0.0; 3]; - - b.spmv_assign_add(&x, &base, &mut y, true); - - assert_eq!(y, base); - } - - #[test] - fn test_from_dense_table_single_element() { - let b = CsrBlock::from_dense_table(&[42.0], 1, 1); - assert_eq!(b.nrows, 1); - assert_eq!(b.ncols, 1); - assert_eq!(b.nnz(), 1); - assert_eq!(b.data[0], 42.0); - assert_eq!(b.indices[0], 0); - } - - #[test] - fn test_transpose_empty() { - let b = CsrBlock::from_dense_table(&[0.0; 6], 2, 3); - let bt = b.transpose(); - assert_eq!(bt.nrows, 3); - assert_eq!(bt.ncols, 2); - assert_eq!(bt.nnz(), 0); - } - - // Chunk *sizing* cannot change the result, so it is deliberately not asserted. - #[test] - fn par_spmv_matches_seq_above_threshold() { - let nrows = PAR_SPMV_THRESHOLD + 2_000; - let ncols = 64; - let half = ncols / 2; - - // Values vary with the row, so a misindexed row is observable. - let mut indptr = vec![0u32; nrows + 1]; - let mut indices = Vec::with_capacity(2 * nrows); - let mut data = Vec::with_capacity(2 * nrows); - for i in 0..nrows { - indices.push((i % half) as u32); - data.push(i as f64 * 0.5 + 1.0); - indices.push((half + i % half) as u32); - data.push(i as f64 * -0.25 - 2.0); - indptr[i + 1] = indptr[i] + 2; - } - let block = CsrBlock { - indptr, - indices, - data, - nrows, - ncols, - }; - - let x: Vec = (0..ncols).map(|j| (j as f64).sin()).collect(); - let base: Vec = (0..nrows).map(|i| i as f64 * 1e-3).collect(); - - let mut y_par = vec![0.0; nrows]; - block.par_spmv_assign_add(&x, &base, &mut y_par); - - let mut y_seq = vec![0.0; nrows]; - block.seq_spmv_assign_add(&x, &base, &mut y_seq); - - assert_eq!(y_par, y_seq); - } -} diff --git a/crates/within/src/domain.rs b/crates/within/src/domain.rs index b9727da5..263fa005 100644 --- a/crates/within/src/domain.rs +++ b/crates/within/src/domain.rs @@ -1,6 +1,5 @@ //! Domain layer: [`Design`] (design-matrix metadata) and factor-pair [`Subdomain`] construction. -pub(crate) mod cross_tab; mod effect; pub(crate) mod factor_pairs; mod level_moments; @@ -8,17 +7,13 @@ mod prepared; mod reparam; mod scatter_split; -pub(crate) use cross_tab::CrossTab; pub(crate) use prepared::{PreparedDesign, PreparedTerm}; use reparam::TermReparam; pub(crate) use scatter_split::{LevelRows, ScatterSplit}; pub use effect::Effect; -pub(crate) use factor_pairs::{ - build_local_domains, CoordinateMap, Grounding, LocalComponent, LocalDomain, MatrixForm, - SddmMatrix, -}; +pub(crate) use factor_pairs::{build_local_domains, LocalDomain}; use crate::channel::Channel; use crate::BuildError; diff --git a/crates/within/src/domain/cross_tab.rs b/crates/within/src/domain/cross_tab.rs deleted file mode 100644 index 458aaa02..00000000 --- a/crates/within/src/domain/cross_tab.rs +++ /dev/null @@ -1,237 +0,0 @@ -//! Cross-tabulation of a channel pair: the bipartite local Gramian. -//! -//! [`CrossTab`] holds `C` as a [`CsrBlock`] plus its transpose, built on first use -//! (rather than assembling the symmetric block matrix), and supports bipartite -//! connected-components splitting and per-component extraction. -//! Levels use the design's compact positions, with a `local_to_global` map into -//! the full coefficient space. - -use std::sync::OnceLock; - -use serde::ser::SerializeStruct; - -use crate::channel::ChannelPair; -use crate::csr_block::{to_u32, CsrBlock}; -use crate::domain::PreparedDesign; - -mod accumulate; -use accumulate::accumulate_cross_block; -pub(crate) use accumulate::PairColumns; - -/// A connected component in a bipartite factor-pair graph, in compact 0-based parent indices. -pub(crate) struct BipartiteComponent { - pub(crate) rows: Vec, - pub(crate) cols: Vec, -} - -/// The off-diagonal block `C` of a channel pair's local Gramian, plus `Cᵀ`. -#[derive(Clone)] -pub(crate) struct CrossTab { - /// CSR(C): row-block rows (n_rows) x col-block cols (n_cols). - pub(crate) c: CsrBlock, - /// CSR(Cᵀ), built on first use so a cover that only reads rows of `C` never pays for it. - ct: OnceLock, -} - -impl CrossTab { - pub(crate) fn new(c: CsrBlock) -> Self { - Self { - c, - ct: OnceLock::new(), - } - } - - pub(crate) fn with_transpose(c: CsrBlock, ct: CsrBlock) -> Self { - debug_assert!(ct.nrows == c.ncols && ct.ncols == c.nrows); - Self { - c, - ct: OnceLock::from(ct), - } - } - - /// Both blocks up front, for a cross-tab the solve will read on every iteration. - pub(crate) fn eager(c: CsrBlock) -> Self { - let ct = c.transpose(); - Self::with_transpose(c, ct) - } - - pub(crate) fn ct(&self) -> &CsrBlock { - self.ct.get_or_init(|| self.c.transpose()) - } - - pub(crate) fn into_parts(self) -> (CsrBlock, CsrBlock) { - self.ct(); - let ct = self.ct.into_inner().expect("transpose was just built"); - (self.c, ct) - } -} - -// The wire carries both blocks, as the derive did, so the format is unchanged. -impl serde::Serialize for CrossTab { - fn serialize(&self, serializer: S) -> Result { - let mut state = serializer.serialize_struct("CrossTab", 2)?; - state.serialize_field("c", &self.c)?; - state.serialize_field("ct", self.ct())?; - state.end() - } -} - -impl<'de> serde::Deserialize<'de> for CrossTab { - fn deserialize>(deserializer: D) -> Result { - #[derive(serde::Deserialize)] - #[serde(rename = "CrossTab")] - struct Wire { - c: CsrBlock, - ct: CsrBlock, - } - // Shapes are the owning solver's decoder's to reject; corrupt bytes must not panic here. - let Wire { c, ct } = Wire::deserialize(deserializer)?; - Ok(Self { - c, - ct: OnceLock::from(ct), - }) - } -} - -impl CrossTab { - /// Number of rows in the row block. - pub(crate) fn n_rows(&self) -> usize { - self.c.nrows - } - - /// Number of rows in the col block. - pub(crate) fn n_cols(&self) -> usize { - self.c.ncols - } - - /// Total number of DOFs (n_rows + n_cols). - pub(crate) fn n_local(&self) -> usize { - self.c.nrows + self.c.ncols - } - - /// Build one channel pair's `C` and its local-to-global DOF map. - pub(crate) fn build_for_pair( - prepared: &PreparedDesign<'_>, - pair: ChannelPair, - ) -> (Self, Vec) { - let design = &prepared.design; - let row_term = &design.terms[pair.rows.term]; - let col_term = &design.terms[pair.cols.term]; - let (n_rows, n_cols) = (row_term.n_levels(), col_term.n_levels()); - - let cross_tab = CrossTab::eager(accumulate_cross_block(prepared, pair, n_rows, n_cols)); - let row_base = row_term.column_dofs(pair.rows.column).start; - let col_base = col_term.column_dofs(pair.cols.column).start; - let local_to_global = (0..n_rows) - .map(|level| to_u32(row_base + level)) - .chain((0..n_cols).map(|level| to_u32(col_base + level))) - .collect(); - - (cross_tab, local_to_global) - } - - /// Symmetric adjacency over local `[q | r]` indexing: q-nodes walk `C`, r-nodes walk `Cᵀ`. - pub(crate) fn neighbors(&self, i: usize) -> impl Iterator + '_ { - let n_rows = self.n_rows(); - let (block, row, off) = if i < n_rows { - (&self.c, i, n_rows) - } else { - (self.ct(), i - n_rows, 0) - }; - block.row(row).map(move |(j, v)| (off + j, v)) - } - - /// Connected components by DFS over [`Self::neighbors`]; O(n_rows + n_cols + nnz). - pub(crate) fn bipartite_connected_components(&self) -> Vec { - let n_rows = self.n_rows(); - let mut visited = vec![false; self.n_local()]; - let mut components = Vec::new(); - let mut stack = Vec::new(); - - for start in 0..self.n_local() { - if visited[start] { - continue; - } - visited[start] = true; - stack.push(start); - let mut rows = Vec::new(); - let mut cols = Vec::new(); - - while let Some(node) = stack.pop() { - if node < n_rows { - rows.push(node); - } else { - cols.push(node - n_rows); - } - for (j, _) in self.neighbors(node) { - if !visited[j] { - visited[j] = true; - stack.push(j); - } - } - } - - rows.sort_unstable(); - cols.sort_unstable(); - components.push(BipartiteComponent { rows, cols }); - } - - components - } - - /// `row_remap`/`col_remap` must arrive all-`u32::MAX` and are reset on exit. - pub(crate) fn extract_component( - &self, - comp: &BipartiteComponent, - row_remap: &mut [u32], - col_remap: &mut [u32], - ) -> Self { - let n_rows = comp.rows.len(); - let n_cols = comp.cols.len(); - debug_assert_eq!(row_remap.len(), self.n_rows()); - debug_assert_eq!(col_remap.len(), self.n_cols()); - - for (new_idx, &old_idx) in comp.rows.iter().enumerate() { - row_remap[old_idx] = to_u32(new_idx); - } - for (new_idx, &old_idx) in comp.cols.iter().enumerate() { - col_remap[old_idx] = to_u32(new_idx); - } - - let mut c_indptr = vec![0u32; n_rows + 1]; - let mut c_indices = Vec::new(); - let mut c_data = Vec::new(); - for (new_row, &old_row) in comp.rows.iter().enumerate() { - let start = self.c.indptr[old_row] as usize; - let end = self.c.indptr[old_row + 1] as usize; - for idx in start..end { - let old_rj = self.c.indices[idx] as usize; - let new_rj = col_remap[old_rj]; - if new_rj != u32::MAX { - c_indices.push(new_rj); - c_data.push(self.c.data[idx]); - } - } - c_indptr[new_row + 1] = to_u32(c_indices.len()); - } - - let c = CsrBlock { - indptr: c_indptr, - indices: c_indices, - data: c_data, - nrows: n_rows, - ncols: n_cols, - }; - for &old_idx in &comp.rows { - row_remap[old_idx] = u32::MAX; - } - for &old_idx in &comp.cols { - col_remap[old_idx] = u32::MAX; - } - - CrossTab::eager(c) - } -} - -#[cfg(test)] -mod tests; diff --git a/crates/within/src/domain/cross_tab/accumulate.rs b/crates/within/src/domain/cross_tab/accumulate.rs deleted file mode 100644 index 2bb042d4..00000000 --- a/crates/within/src/domain/cross_tab/accumulate.rs +++ /dev/null @@ -1,270 +0,0 @@ -//! Observation accumulation kernels for [`CrossTab`](super::CrossTab) construction. -//! -//! Both the dense and sparse paths scan [`PairColumns::observations`] once, adding -//! `w·l_row·l_col` to each cell, where `l` is the channel's loading, so slope channels -//! yield signed cells. -//! Paths are generic over [`Loading`] and monomorphized per pair: intercept -//! channels pass [`Unit`], whose `l ≡ 1` folds the loading math away, so plain -//! pairs keep the pre-slope codegen. - -use crate::channel::ChannelPair; -use crate::csr_block::CsrBlock; -use crate::domain::{row_weight, PreparedDesign}; -use crate::linalg::add_compensated; - -use super::to_u32; - -/// Prefer sparse above this accumulator size when it reduces transient memory. -const DENSE_TABLE_MAX_BYTES: usize = 40_000_000; - -/// A channel's per-observation loading. -pub(crate) trait Loading: Copy { - fn at(self, uid: usize) -> f64; -} - -/// Intercept loading `l ≡ 1`; LLVM folds the resulting `x · 1.0` away. -#[derive(Clone, Copy)] -pub(super) struct Unit; - -impl Loading for Unit { - #[inline] - fn at(self, _uid: usize) -> f64 { - 1.0 - } -} - -impl Loading for &[f64] { - #[inline] - fn at(self, uid: usize) -> f64 { - self[uid] - } -} - -/// `None` is an intercept's `l ≡ 1`, for a reader that needs no specialized kernel. -impl Loading for Option<&[f64]> { - #[inline] - fn at(self, uid: usize) -> f64 { - self.map_or(1.0, |z| z[uid]) - } -} - -/// One weighted observation of a channel pair, at its cell with both loadings. -pub(crate) struct Observation { - pub(crate) row: usize, - pub(crate) col: usize, - pub(crate) x: f64, - pub(crate) y: f64, - pub(crate) w: f64, -} - -/// Per-observation input columns backing one channel pair: level codes, loadings, and weights. -#[derive(Clone, Copy)] -pub(crate) struct PairColumns<'a, Lq: Loading, Lr: Loading> { - pub(super) row_levels: &'a [u32], - pub(super) col_levels: &'a [u32], - pub(super) row_load: Lq, - pub(super) col_load: Lr, - pub(super) sqrt_weights: Option<&'a [f64]>, -} - -impl<'a> PairColumns<'a, Option<&'a [f64]>, Option<&'a [f64]>> { - pub(crate) fn new(prepared: &'a PreparedDesign<'_>, pair: ChannelPair) -> Self { - let (rows, cols) = (prepared.term(pair.rows.term), prepared.term(pair.cols.term)); - Self { - row_levels: rows.term.levels(), - col_levels: cols.term.levels(), - row_load: rows.loading(pair.rows.column), - col_load: cols.loading(pair.cols.column), - sqrt_weights: prepared.sqrt_weights(), - } - } -} - -impl PairColumns<'_, Lq, Lr> { - #[inline] - pub(crate) fn observations(&self) -> impl Iterator + '_ { - (0..self.row_levels.len()).map(move |uid| Observation { - row: self.row_levels[uid] as usize, - col: self.col_levels[uid] as usize, - x: self.row_load.at(uid), - y: self.col_load.at(uid), - w: row_weight(self.sqrt_weights, uid), - }) - } -} - -/// Accumulate into `C`, dispatching dense or sparse by peak transient memory. -pub(super) fn accumulate_cross_block( - prepared: &PreparedDesign<'_>, - pair: ChannelPair, - n_rows: usize, - n_cols: usize, -) -> CsrBlock { - let design = &prepared.design; - let sqrt_weights = prepared.sqrt_weights(); - // Dispatching on cell count alone would pick sparse where it uses MORE memory. - let table_size = n_rows.saturating_mul(n_cols); - let dense_cost = table_size.saturating_mul(16); - let sparse_cost = design.n_obs.saturating_mul(12); - let go_sparse = dense_cost > DENSE_TABLE_MAX_BYTES && sparse_cost < dense_cost; - - let (rows, cols) = (prepared.term(pair.rows.term), prepared.term(pair.cols.term)); - let (row_levels, col_levels) = (rows.term.levels(), cols.term.levels()); - // One arm per loading combination; closures aren't generic, so the literals repeat. - match ( - rows.loading(pair.rows.column), - cols.loading(pair.cols.column), - ) { - (None, None) => accumulate( - PairColumns { - row_levels, - col_levels, - row_load: Unit, - col_load: Unit, - sqrt_weights, - }, - n_rows, - n_cols, - go_sparse, - ), - (Some(zq), None) => accumulate( - PairColumns { - row_levels, - col_levels, - row_load: zq, - col_load: Unit, - sqrt_weights, - }, - n_rows, - n_cols, - go_sparse, - ), - (None, Some(zr)) => accumulate( - PairColumns { - row_levels, - col_levels, - row_load: Unit, - col_load: zr, - sqrt_weights, - }, - n_rows, - n_cols, - go_sparse, - ), - (Some(zq), Some(zr)) => accumulate( - PairColumns { - row_levels, - col_levels, - row_load: zq, - col_load: zr, - sqrt_weights, - }, - n_rows, - n_cols, - go_sparse, - ), - } -} - -/// Size-dispatched accumulation for one monomorphized loading combination. -fn accumulate( - cols: PairColumns<'_, Lq, Lr>, - n_rows: usize, - n_cols: usize, - go_sparse: bool, -) -> CsrBlock { - if go_sparse { - accumulate_sparse_cross_block(cols, n_rows, n_cols) - } else { - accumulate_dense_cross_block(cols, n_rows, n_cols) - } -} - -/// Dense path: flat `n_rows * n_cols` table with O(1) accumulation per observation. -pub(super) fn accumulate_dense_cross_block( - cols: PairColumns<'_, Lq, Lr>, - n_rows: usize, - n_cols: usize, -) -> CsrBlock { - let mut table = vec![0.0f64; n_rows * n_cols]; - let mut corrections = vec![0.0f64; table.len()]; - - for o in cols.observations() { - debug_assert!(o.row < n_rows && o.col < n_cols); - let cell = o.row * n_cols + o.col; - add_compensated(&mut table[cell], &mut corrections[cell], o.w * o.x * o.y); - } - - for (cell, correction) in table.iter_mut().zip(corrections) { - *cell += correction; - } - CsrBlock::from_dense_table(&table, n_rows, n_cols) -} - -/// Sparse path: bucket by row, then dedup each row through a dense `n_cols` workspace. -pub(super) fn accumulate_sparse_cross_block( - cols: PairColumns<'_, Lq, Lr>, - n_rows: usize, - n_cols: usize, -) -> CsrBlock { - let mut row_counts = vec![0u32; n_rows]; - for &row in cols.row_levels { - row_counts[row as usize] += 1; - } - - let mut bucket_indptr = vec![0u32; n_rows + 1]; - for i in 0..n_rows { - bucket_indptr[i + 1] = bucket_indptr[i] + row_counts[i]; - } - let total_entries = bucket_indptr[n_rows] as usize; - - let mut bucket_cols = vec![0u32; total_entries]; - let mut bucket_vals = vec![0.0f64; total_entries]; - let mut cursor = bucket_indptr[..n_rows].to_vec(); - for o in cols.observations() { - let pos = cursor[o.row] as usize; - bucket_cols[pos] = to_u32(o.col); - bucket_vals[pos] = o.w * o.x * o.y; - cursor[o.row] += 1; - } - - // A signed cell cancelling to 0.0 mid-row re-pushes its column; the duplicate is harmless. - let mut work = vec![0.0f64; n_cols]; - let mut corrections = vec![0.0f64; n_cols]; - let mut touched: Vec = Vec::new(); - let mut c_indptr = vec![0u32; n_rows + 1]; - let mut c_indices = Vec::new(); - let mut c_data = Vec::new(); - - for row in 0..n_rows { - let start = bucket_indptr[row] as usize; - let end = bucket_indptr[row + 1] as usize; - for idx in start..end { - let col = bucket_cols[idx] as usize; - if work[col] == 0.0 && corrections[col] == 0.0 { - touched.push(to_u32(col)); - } - add_compensated(&mut work[col], &mut corrections[col], bucket_vals[idx]); - } - touched.sort_unstable(); - for &col in &touched { - let v = work[col as usize] + corrections[col as usize]; - if v != 0.0 { - c_indices.push(col); - c_data.push(v); - } - work[col as usize] = 0.0; - corrections[col as usize] = 0.0; - } - c_indptr[row + 1] = to_u32(c_indices.len()); - touched.clear(); - } - - CsrBlock { - indptr: c_indptr, - indices: c_indices, - data: c_data, - nrows: n_rows, - ncols: n_cols, - } -} diff --git a/crates/within/src/domain/cross_tab/tests.rs b/crates/within/src/domain/cross_tab/tests.rs deleted file mode 100644 index 8990e8cf..00000000 --- a/crates/within/src/domain/cross_tab/tests.rs +++ /dev/null @@ -1,277 +0,0 @@ -use proptest::prelude::*; - -use super::accumulate::{ - accumulate_dense_cross_block, accumulate_sparse_cross_block, PairColumns, Unit, -}; -use super::CrossTab; -use crate::channel::{Channel, ChannelPair}; -use crate::csr_block::CsrBlock; -use crate::domain::{Design, Effect, PreparedDesign}; - -impl CrossTab { - pub(crate) fn from_dense_for_test(table: &[f64], n_rows: usize, n_cols: usize) -> Self { - let c = CsrBlock::from_dense_table(table, n_rows, n_cols); - Self::eager(c) - } -} - -/// Terms 0 and 1 paired on their intercept channels (plain cross-tab). -const INTERCEPT_PAIR: ChannelPair = ChannelPair { - rows: Channel { term: 0, column: 0 }, - cols: Channel { term: 1, column: 0 }, -}; - -fn design_of(columns: Vec>) -> PreparedDesign<'static> { - PreparedDesign::from_levels_for_test(columns) -} - -#[test] -fn over_threshold_pair_matches_the_dense_kernel() { - // 2,500 × 2,500 levels exceed the dense cap; unsorted rows see each column twice, descending. - let n_levels = 2500usize; - let fa: Vec = (0..8 * n_levels).map(|i| (i % n_levels) as u32).collect(); - let fb: Vec = (0..8 * n_levels) - .map(|i| ((i % n_levels * 7 + (3 - i / n_levels % 4) * 13) % n_levels) as u32) - .collect(); - let effects = [&fa, &fb].map(|levels| Effect::new(levels, true, []).unwrap()); - let design = PreparedDesign::unweighted_for_test(Design::new_unsorted(effects).unwrap()); - let (ct, _) = CrossTab::build_for_pair(&design, INTERCEPT_PAIR); - - let terms = &design.design.terms; - let cols = PairColumns { - row_levels: terms[0].levels(), - col_levels: terms[1].levels(), - row_load: Unit, - col_load: Unit, - sqrt_weights: None, - }; - let dense = accumulate_dense_cross_block(cols, terms[0].n_levels(), terms[1].n_levels()); - assert_eq!(ct.c.indptr, dense.indptr); - assert_eq!(ct.c.indices, dense.indices); - assert_eq!(ct.c.data, dense.data); -} - -#[test] -fn test_extract_component_two_components() { - // Two disconnected bipartite components: q/r levels {0,1} and {2,3}. - let fa = vec![0u32, 0, 1, 1, 2, 2, 3, 3]; - let fb = vec![0u32, 1, 0, 1, 2, 3, 2, 3]; - let design = design_of(vec![fa, fb]); - let (ct, _) = CrossTab::build_for_pair(&design, INTERCEPT_PAIR); - - let components = ct.bipartite_connected_components(); - assert_eq!(components.len(), 2, "should have 2 connected components"); - - // Reusable remap buffers, reset by `extract_component` between components. - let mut row_remap = vec![u32::MAX; ct.n_rows()]; - let mut col_remap = vec![u32::MAX; ct.n_cols()]; - - // Sort components by their smallest q-index for deterministic comparison. - let mut comps: Vec<_> = components.iter().collect(); - comps.sort_by_key(|c| c.rows[0]); - - let comp_a = comps[0]; - let comp_b = comps[1]; - - assert_eq!(comp_a.rows, vec![0, 1], "component A row indices"); - assert_eq!(comp_a.cols, vec![0, 1], "component A col indices"); - assert_eq!(comp_b.rows, vec![2, 3], "component B row indices"); - assert_eq!(comp_b.cols, vec![2, 3], "component B col indices"); - - // Extract component A and verify its sub-CrossTab. - let sub_a = ct.extract_component(comp_a, &mut row_remap, &mut col_remap); - assert_eq!(sub_a.n_rows(), 2, "component A: n_rows=2"); - assert_eq!(sub_a.n_cols(), 2, "component A: n_cols=2"); - - // Column indices in sub_a.c should be 0-based (0..n_cols for component A = 0..2). - let max_col_a = sub_a.c.indices.iter().copied().max().unwrap_or(0); - assert!( - (max_col_a as usize) < sub_a.n_cols(), - "sub_a C column indices should be 0-based < n_cols={}", - sub_a.n_cols() - ); - - // C^T of sub_a should equal the exact transpose of sub_a.c. - let ct_t = sub_a.c.transpose(); - assert_eq!( - ct_t.indptr, - sub_a.ct().indptr, - "sub_a: ct.indptr should equal transpose(c).indptr" - ); - assert_eq!( - ct_t.indices, - sub_a.ct().indices, - "sub_a: ct.indices should equal transpose(c).indices" - ); - for (a, b) in ct_t.data.iter().zip(&sub_a.ct().data) { - assert!( - (a - b).abs() < 1e-12, - "sub_a: ct.data should equal transpose(c).data" - ); - } - - // Extract component B and verify its sub-CrossTab. - let sub_b = ct.extract_component(comp_b, &mut row_remap, &mut col_remap); - assert_eq!(sub_b.n_rows(), 2, "component B: n_rows=2"); - assert_eq!(sub_b.n_cols(), 2, "component B: n_cols=2"); - - // Column indices in sub_b.c should be 0-based. - let max_col_b = sub_b.c.indices.iter().copied().max().unwrap_or(0); - assert!( - (max_col_b as usize) < sub_b.n_cols(), - "sub_b C column indices should be 0-based < n_cols={}", - sub_b.n_cols() - ); - - // The two components should have the same structure (symmetric design). - assert_eq!( - sub_a.c.indptr, sub_b.c.indptr, - "symmetric design: sub_a and sub_b should have same C structure" - ); -} - -proptest! { - #![proptest_config(ProptestConfig::with_cases(10))] - - #[test] - fn prop_bipartite_components_partition( - n_rows in 2usize..=8, - n_cols in 2usize..=8, - n_obs in 4usize..=30, - seed in 0u64..1000, - ) { - // Generate observations using a deterministic pseudo-random pattern. - let mut fa: Vec = Vec::with_capacity(n_obs); - let mut fb: Vec = Vec::with_capacity(n_obs); - let mut s = seed; - for _ in 0..n_obs { - // LCG: x_{n+1} = (a * x_n + c) mod m - s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); - fa.push((s % n_rows as u64) as u32); - s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); - fb.push((s % n_cols as u64) as u32); - } - - let design = design_of(vec![fa, fb]); - let (ct, _) = CrossTab::build_for_pair(&design, INTERCEPT_PAIR); - - let components = ct.bipartite_connected_components(); - - // Collect all row indices and col indices across components. - let mut all_rows: Vec = components.iter().flat_map(|c| c.rows.iter().copied()).collect(); - let mut all_cols: Vec = components.iter().flat_map(|c| c.cols.iter().copied()).collect(); - all_rows.sort_unstable(); - all_cols.sort_unstable(); - - // Union should cover every compact level position. - let expected_rows: Vec = (0..ct.n_rows()).collect(); - let expected_cols: Vec = (0..ct.n_cols()).collect(); - prop_assert_eq!(&all_rows, &expected_rows, "row indices should cover 0..n_rows={}", ct.n_rows()); - prop_assert_eq!(&all_cols, &expected_cols, "col indices should cover 0..n_cols={}", ct.n_cols()); - - // Indices within each component should be sorted. - for (ci, comp) in components.iter().enumerate() { - prop_assert!( - comp.rows.windows(2).all(|w| w[0] < w[1]), - "component {ci}: rows should be sorted" - ); - prop_assert!( - comp.cols.windows(2).all(|w| w[0] < w[1]), - "component {ci}: cols should be sorted" - ); - } - - // Index sets should be disjoint between components. - let mut rows_seen = std::collections::HashSet::new(); - let mut r_seen = std::collections::HashSet::new(); - for (ci, comp) in components.iter().enumerate() { - for &qi in &comp.rows { - prop_assert!(rows_seen.insert(qi), "component {ci}: q-index {qi} appears in multiple components"); - } - for &ri in &comp.cols { - prop_assert!(r_seen.insert(ri), "component {ci}: r-index {ri} appears in multiple components"); - } - } - } -} - -#[test] -fn design_contains_every_compact_level_position() { - // Cover identity, bounded-gappy, and wide-sparse encodings. - let design = design_of(vec![ - vec![0u32, 1, 2, 0, 1, 2], - vec![0u32, 2, 4, 0, 2, 4], - vec![7u32, u32::MAX, 7, u32::MAX, 7, u32::MAX], - ]); - - for (term, meta) in design.design.terms.iter().enumerate() { - let mut observed = vec![false; meta.n_levels()]; - for &level in meta.levels() { - observed[level as usize] = true; - } - assert!( - observed.into_iter().all(|is_observed| is_observed), - "term {term} contains an unobserved internal position" - ); - } -} - -#[test] -fn dense_and_sparse_paths_agree_on_signed_data() { - // Cell (f=0,g=0) crosses 0.0 mid-row and (f=0,g=1) cancels to 0.0; both paths must drop it. - let f = [0u32, 0, 0, 0, 0, 1]; - let z = [1.0, -1.0, 2.0, 3.0, -3.0, 4.0]; - let g = [0u32, 0, 0, 1, 1, 0]; - let effects = vec![ - Effect::new(&f, true, [&z[..]]).unwrap(), - Effect::new(&g, true, []).unwrap(), - ]; - let design = Design::new(effects).unwrap(); - let pair = ChannelPair { - rows: Channel { term: 0, column: 1 }, - cols: Channel { term: 1, column: 0 }, - }; - let cols = PairColumns { - row_levels: design.terms[0].levels(), - col_levels: design.terms[1].levels(), - row_load: design.terms[0].raw_slopes().next().unwrap(), - col_load: Unit, - sqrt_weights: None, - }; - let n_rows = design.terms[pair.rows.term].n_levels(); - let n_cols = design.terms[pair.cols.term].n_levels(); - let c_dense = accumulate_dense_cross_block(cols, n_rows, n_cols); - let c_sparse = accumulate_sparse_cross_block(cols, n_rows, n_cols); - - // Bit-exact parity: identical per-cell addition order in both paths. - assert_eq!(c_dense.indptr, c_sparse.indptr); - assert_eq!(c_dense.indices, c_sparse.indices); - assert_eq!(c_dense.data, c_sparse.data); - - // Row f=0 keeps only cell (0,0) = 2.0; the exact-0.0 cell (0,1) is gone. - assert_eq!(&c_dense.indptr, &[0, 1, 2]); - assert_eq!(c_dense.indices[0], 0); - assert_eq!(c_dense.data[0], 2.0); -} - -#[test] -fn both_accumulators_recover_signed_cells_and_reset_compensation() { - let rows = [0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1]; - let cols = [0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0]; - let loading = [1e16, 1., -1e16, 1., -1., 0., 2., 0., 1., 1e16, -1e16]; - let pair = PairColumns { - row_levels: &rows, - col_levels: &cols, - row_load: &loading[..], - col_load: Unit, - sqrt_weights: None, - }; - for actual in [ - accumulate_dense_cross_block(pair, 2, 2), - accumulate_sparse_cross_block(pair, 2, 2), - ] { - assert_eq!(actual.indptr, [0, 2, 3]); - assert_eq!(actual.indices, [0, 1, 0]); - assert_eq!(actual.data, [1., 2., 1.]); - } -} diff --git a/crates/within/src/domain/factor_pairs.rs b/crates/within/src/domain/factor_pairs.rs index b96e0f80..ec55e8cd 100644 --- a/crates/within/src/domain/factor_pairs.rs +++ b/crates/within/src/domain/factor_pairs.rs @@ -1,43 +1,32 @@ //! Channel-pair subdomain construction. //! -//! Each cross-factor channel pair becomes a Schwarz subdomain (one per -//! connected component of its bipartite cross-tab; isolated levels are batched -//! into diagonal subdomains per side). Overlap is handled by +//! Each cross-factor channel pair hands approx-chol its Gram factor `V`, one column per +//! observation, and each irreducible block of the factor becomes a Schwarz subdomain +//! (isolated levels are batched into diagonal subdomains per side). Overlap is handled by //! partition-of-unity weights — see [`schwarz_precond::domain`] for the math. //! //! Entry point: [`build_local_domains`]. +use approx_chol::input::{BipartiteFactorWidth2, FactorWidth1}; +use approx_chol::{BlockFactor, ReduceError}; use schwarz_precond::{PartitionWeights, SubdomainCore}; +use crate::block_solver::BlockSolver; use crate::channel::{Channel, ChannelPair}; -use crate::config::LocalSolverConfig; -use crate::csr_block::to_u32; +use crate::config::{LocalSolverConfig, ScalingFailure}; +use crate::domain::Column; use crate::{BuildError, BuildWarning}; -use super::{CrossTab, PreparedDesign}; - -mod grounding; -mod sddm; -use crate::domain::cross_tab::{BipartiteComponent, PairColumns}; -use crate::domain::Column; -use grounding::scaled_groundings; -use sddm::{convert, NotScalable}; -pub(crate) use sddm::{CoordinateMap, Grounding, LocalComponent, MatrixForm, SddmMatrix}; +use super::{PreparedDesign, PreparedTerm}; /// Keeps a batch below the local solver's inner-parallel threshold and typical scratch size. const ISOLATED_LEVEL_CHUNK: usize = 4096; -#[derive(Clone, Copy, Debug, PartialEq, Eq)] -enum ComponentClass { - KnownLaplacian, - General, -} - -/// A factor-pair Schwarz domain paired with its validated local operator. +/// A factor-pair Schwarz domain paired with its factored local operator. #[derive(Clone)] pub(crate) struct LocalDomain { pub(crate) core: SubdomainCore, - pub(crate) component: LocalComponent, + pub(crate) solver: BlockSolver, } /// Same-factor channel pairs are exactly orthogonal after whitening, so never enumerated. @@ -66,17 +55,7 @@ pub(crate) fn build_local_domains( .collect(); let per_pair: Vec<(Vec, Vec)> = pairs .par_iter() - .map(|&pair| { - let (full_ct, l2g) = CrossTab::build_for_pair(prepared, pair); - let class = if design.column(pair.rows) == Column::Intercept - && design.column(pair.cols) == Column::Intercept - { - ComponentClass::KnownLaplacian - } else { - ComponentClass::General - }; - split_into_subdomains(prepared, pair, class, full_ct, &l2g, config) - }) + .map(|&pair| factor_pair(prepared, pair, config)) .collect::>()?; let mut domain_pairs = Vec::new(); let mut warnings = Vec::new(); @@ -96,119 +75,91 @@ pub(crate) fn build_local_domains( Ok((domain_pairs, warnings)) } -/// Dead singletons (zero diagonal, an exact-zero design column) produce no subdomain. -fn split_into_subdomains( +/// `G = V Vᵀ`, observation `k` adding `√w_k (z_q(k) e_{a_k} + z_r(k) e_{n_q + b_k})`. +fn factor_pair( prepared: &PreparedDesign<'_>, pair: ChannelPair, - class: ComponentClass, - full_ct: CrossTab, - l2g: &[u32], config: &LocalSolverConfig, ) -> Result<(Vec, Vec), BuildError> { - let row_diag = prepared.channel_diagonal(pair.rows); - let col_diag = prepared.channel_diagonal(pair.cols); - debug_assert_eq!( - (row_diag.len(), col_diag.len()), - (full_ct.n_rows(), full_ct.n_cols()) - ); - let n_rows_full = full_ct.n_rows(); - let mut components = Vec::new(); - let (mut isolated_rows, mut isolated_cols) = (Vec::new(), Vec::new()); - for comp in full_ct.bipartite_connected_components() { - match (&comp.rows[..], &comp.cols[..]) { - (&[i], []) if row_diag[i] != 0.0 => isolated_rows.push(i), - ([], &[j]) if col_diag[j] != 0.0 => isolated_cols.push(j), - _ => components.push(comp), - } + let (rows, cols) = (prepared.term(pair.rows.term), prepared.term(pair.cols.term)); + fn side<'p>(term: &PreparedTerm<'p>, column: usize) -> FactorWidth1<'p, f64> { + FactorWidth1::new( + term.term.n_levels(), + term.term.levels(), + term.loading(column), + ) + .expect("levels index their term") } - // Isolated levels have a diagonal Gram, so batching changes no local solve, only overhead. - components.extend( - isolated_rows - .chunks(ISOLATED_LEVEL_CHUNK) - .map(|rows| BipartiteComponent { - rows: rows.to_vec(), - cols: Vec::new(), - }), - ); - components.extend( - isolated_cols - .chunks(ISOLATED_LEVEL_CHUNK) - .map(|cols| BipartiteComponent { - rows: Vec::new(), - cols: cols.to_vec(), - }), - ); - - let cross_tabs: Vec = if components.len() == 1 { - vec![full_ct] - } else { - let mut row_remap = vec![u32::MAX; full_ct.n_rows()]; - let mut col_remap = vec![u32::MAX; full_ct.n_cols()]; - components - .iter() - .map(|comp| full_ct.extract_component(comp, &mut row_remap, &mut col_remap)) - .collect() - }; + let (left, right) = (side(&rows, pair.rows.column), side(&cols, pair.cols.column)); + let v = match prepared.sqrt_weights() { + Some(s) => BipartiteFactorWidth2::scaled(left, right, s), + None => BipartiteFactorWidth2::new(left, right), + } + .expect("both sides read the same observations"); + let factor = approx_chol::factorize_with(v, config.factor_config(), config.scaling_policy()?) + .map_err(|e| match e { + approx_chol::Error::Reduce(ReduceError::Unscalable { .. }) => { + BuildError::UnscalableComponent { pair } + } + e => BuildError::LocalSolverBuild(e.to_string()), + })?; - // Only folded components are tested: a frustrated block is nonsingular, so it never floats. - let mut levels: Option>> = None; - let mut domains = Vec::with_capacity(components.len()); let mut warnings = Vec::new(); - for (comp, comp_ct) in components.iter().zip(cross_tabs) { - let comp_diag: Vec = comp - .rows - .iter() - .map(|&i| row_diag[i]) - .chain(comp.cols.iter().map(|&i| col_diag[i])) - .collect(); - if comp_diag.iter().all(|&v| v == 0.0) { - continue; - } - // Pair-local `[rows | cols]` positions, mapped to globals once oriented. - let comp_indices: Vec = comp - .rows - .iter() - .map(|&i| to_u32(i)) - .chain(comp.cols.iter().map(|&j| to_u32(n_rows_full + j))) - .collect(); - let (comp_ct, comp_diag, mut comp_indices) = - sddm::orient_for_elimination(comp_ct, comp_diag, comp_indices); - let (component, uncertified) = convert(comp_ct, comp_diag, class, &config.scaling) - .map_err(|NotScalable| BuildError::UnscalableComponent { pair })?; - if let Some(uncertified) = uncertified { - warnings.push(BuildWarning::UnscalableComponent { + for report in factor + .blocks() + .iter() + .filter_map(BlockFactor::scaling_report) + .filter(|report| report.violation > config.scaling.tolerance) + { + match config.scaling.on_failure { + ScalingFailure::Warn => warnings.push(BuildWarning::UnscalableComponent { pair, - iterations: uncertified.iterations, - violation: uncertified.violation, - }); - } - if class == ComponentClass::General && component.form == MatrixForm::Laplacian { - let levels = levels.get_or_insert_with(|| vec![None; n_rows_full + col_diag.len()]); - let mut factors = vec![1.0; comp_indices.len()]; - component - .coordinates - .fold(&mut factors, component.matrix.n_eliminated()); - for (&p, &f) in comp_indices.iter().zip(&factors) { - levels[p as usize] = Some((domains.len(), f)); - } - } - for index in &mut comp_indices { - *index = l2g[*index as usize]; + iterations: report.iterations, + violation: report.violation, + }), + ScalingFailure::Error => return Err(BuildError::UnscalableComponent { pair }), } - domains.push(LocalDomain { - core: schwarz_precond::SubdomainCore::uniform(comp_indices), - component, - }); } - // A folded component's dropped surplus is `fᵀGf` exactly, so only those read the observations. - if let Some(levels) = levels { - let columns = PairColumns::new(prepared, pair); - let groundings = scaled_groundings(&columns, &levels, n_rows_full, domains.len()); - for (ld, grounding) in domains.iter_mut().zip(groundings) { - if grounding == Grounding::Floating { - ld.component.matrix.float(); + let n_rows = rows.term.n_levels(); + let row_base = rows.term.column_dofs(pair.rows.column).start; + let col_base = cols.term.column_dofs(pair.cols.column).start; + let global = |vertex: usize| { + let dof = if vertex < n_rows { + row_base + vertex + } else { + col_base + vertex - n_rows + }; + u32::try_from(dof).expect("the DOF space fits u32") + }; + let domain = |blocks: Vec| LocalDomain { + core: SubdomainCore::uniform( + blocks + .iter() + .flat_map(BlockFactor::vertices) + .map(global) + .collect(), + ), + solver: BlockSolver::new(blocks), + }; + + let mut domains = Vec::new(); + let mut isolated = [Vec::new(), Vec::new()]; + for block in factor.into_blocks() { + match block.vertices().len() { + 1 => isolated[usize::from(block.vertices().all(|v| v >= n_rows))].push(block), + _ => domains.push(domain(vec![block])), + } + } + // Isolated levels have a diagonal Gram, so batching changes no local solve, only overhead. + for side in isolated { + let mut side = side.into_iter(); + loop { + let batch: Vec = side.by_ref().take(ISOLATED_LEVEL_CHUNK).collect(); + if batch.is_empty() { + break; } + domains.push(domain(batch)); } } Ok((domains, warnings)) @@ -341,37 +292,6 @@ mod tests { ); } - /// `z = a_row · b_col` aliases the intercept; a floated block keeps no surplus. - #[rstest::rstest] - #[case::aliased(0.0, Grounding::Floating)] - #[case::identified(1e-4, Grounding::Grounded)] - fn a_folded_alias_floats_without_surplus(#[case] delta: f64, #[case] expected: Grounding) { - let (rows, cols) = ([0u32, 0, 1, 1, 2, 2], [0u32, 1, 0, 1, 0, 1]); - let z = [1.0, 3.0, -2.0, -6.0 * (1.0 + delta), 0.5, 1.5]; - let weights = [0.3, 2.0, 1.0, 7.5, 0.01, 4.0]; - let design = Design::new(vec![ - Effect::new(&rows, false, [&z[..]]).expect("slope effect"), - Effect::new(&cols, true, []).expect("plain effect"), - ]) - .expect("valid design"); - let prepared = PreparedDesign::new(design, Some(&weights)).expect("valid weights"); - let (domains, _) = - build_local_domains(&prepared, &LocalSolverConfig::default()).expect("domains"); - let [ld] = &domains[..] else { - panic!("one component expected, got {}", domains.len()) - }; - let matrix = &ld.component.matrix; - assert_eq!(ld.component.form, MatrixForm::Laplacian); - assert_eq!(matrix.grounding, expected); - if expected == Grounding::Floating { - assert!(matrix.ground_edges.iter().all(|&g| g == 0.0)); - for (i, &d) in matrix.diagonal.iter().enumerate() { - let sum: f64 = matrix.cross_tab.neighbors(i).map(|(_, v)| v.abs()).sum(); - assert_eq!(d, sum, "row {i}"); - } - } - } - #[test] fn isolated_levels_batch_per_side_and_dead_levels_drop() { // ±z within one firm cancels a cross cell to exactly zero; worker 4's zero z is dead. diff --git a/crates/within/src/domain/factor_pairs/grounding.rs b/crates/within/src/domain/factor_pairs/grounding.rs deleted file mode 100644 index 7de977d5..00000000 --- a/crates/within/src/domain/factor_pairs/grounding.rs +++ /dev/null @@ -1,53 +0,0 @@ -//! A folded slope component's [`Grounding`], read from its observations, not its Gram. - -use crate::domain::cross_tab::PairColumns; - -use super::Grounding; - -// 10u: a scaling that misses the kernel leaves surplus that floating would drop. -const SCALED_FLOATING_BOUND: f64 = 5.0 * f64::EPSILON; - -/// `Floating` iff a tested component's scaling `f` has `fᵀGf ≤ 10u·fᵀDf`, per observation. -pub(super) fn scaled_groundings( - columns: &PairColumns<'_, Option<&[f64]>, Option<&[f64]>>, - // Per level `[rows | cols]` of the pair: the tested component and its factor. - levels: &[Option<(usize, f64)>], - n_rows: usize, - n_components: usize, -) -> Vec { - let mut residual = vec![0.0; n_components]; - let mut diagonal = vec![0.0; n_components]; - for o in columns.observations().filter(|o| o.w != 0.0) { - match (levels[o.row], levels[n_rows + o.col]) { - (Some((k, fi)), Some((l, fj))) if k == l => { - let (p, q) = (fi * o.x, fj * o.y); - residual[k] += o.w * (p + q) * (p + q); - diagonal[k] += o.w * (p * p + q * q); - } - (row, col) => { - for (k, p) in [ - row.map(|(k, f)| (k, f * o.x)), - col.map(|(k, f)| (k, f * o.y)), - ] - .into_iter() - .flatten() - { - residual[k] += o.w * p * p; - diagonal[k] += o.w * p * p; - } - } - } - } - residual - .iter() - .zip(&diagonal) - .map(|(&r, &d)| { - // An untested component reads `d = 0`: it never floats. - if d > 0.0 && d.is_finite() && r <= SCALED_FLOATING_BOUND * d { - Grounding::Floating - } else { - Grounding::Grounded - } - }) - .collect() -} diff --git a/crates/within/src/domain/factor_pairs/sddm.rs b/crates/within/src/domain/factor_pairs/sddm.rs deleted file mode 100644 index 60e33aaf..00000000 --- a/crates/within/src/domain/factor_pairs/sddm.rs +++ /dev/null @@ -1,383 +0,0 @@ -//! Validated conversion from a signed bipartite Gram block to SDDM form. -//! -//! One signed diagonal coordinate map `S` makes `S A S` a Z-matrix and scales -//! it to diagonal dominance. Surplus is retained as explicit ground-edge -//! weight, so downstream reduction never has to infer rank or grounding from -//! rounded diagonals. Plain components take a fast path: their Laplacian claim -//! is validated and adopted without sign folding, scaling, or a data rewrite. -//! Frustrated components — where no signature exists — keep their single -//! *signed* matrix (dominance-scaled but not folded to a Z-matrix) and carry -//! a [`MatrixForm::SignedPendingCover`] marker; the Gremban cover that makes the -//! reduced Schur SDDM is built transiently at factor time (see -//! [`crate::block_elim`]), so the stored matrix stays single-sized. - -use super::ComponentClass; -use crate::config::{ScalingConfig, ScalingFailure}; -use crate::domain::CrossTab; - -mod scaling; -use scaling::dominance_scaling; - -const LAPLACIAN_VALIDATION_ULPS: f64 = 64.0; - -/// Gauge of a reduced system: `Floating` anchors one node, `Grounded` factors the full complement. -#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, serde::Serialize, serde::Deserialize)] -pub(crate) enum Grounding { - // Keep this order: postcard encodes by declaration order and the fixture pins Floating = 0. - #[default] - Floating, - Grounded, -} - -/// Which form a component's signature left its matrix in. -#[derive(Clone, Copy, Debug, PartialEq, Eq)] -pub(crate) enum MatrixForm { - /// Folded to a Z-matrix: the minor reduces directly. - Laplacian, - /// Kept signed; a Gremban cover is built at factor time so the stored matrix never doubles. - SignedPendingCover, -} - -/// Map between original and SDDM coordinates; `Canonical` negates the col block, not identity. -#[derive(Clone, Debug, Default, serde::Serialize, serde::Deserialize)] -pub(crate) enum CoordinateMap { - // Keep this order: postcard encodes by declaration order and the fixture pins Canonical = 0. - #[default] - Canonical, - /// Diagonal congruence factors (sign · scale), leaving a frustrated component signed. - Scaled(Box<[f64]>), -} - -impl CoordinateMap { - /// `values` spans the internal system, split at `n_eliminated`. - pub(crate) fn fold(&self, values: &mut [f64], n_eliminated: usize) { - match self { - CoordinateMap::Canonical => { - for value in &mut values[n_eliminated..] { - *value = -*value; - } - } - CoordinateMap::Scaled(factors) => { - debug_assert_eq!(values.len(), factors.len()); - for (value, &factor) in values.iter_mut().zip(factors.iter()) { - *value *= factor; - } - } - } - } - - /// The diagonal maps are involutions up to scale, so they re-apply `fold`. - pub(crate) fn unfold(&self, values: &mut [f64], n_eliminated: usize) { - self.fold(values, n_eliminated); - } -} - -/// Orient so the larger block is eliminated; arrays arrive `[rows | cols]`, leave eliminated-major. -pub(crate) fn orient_for_elimination( - cross_tab: CrossTab, - mut diagonal: Vec, - mut indices: Vec, -) -> (CrossTab, Vec, Vec) { - debug_assert_eq!(indices.len(), cross_tab.n_local()); - debug_assert_eq!(diagonal.len(), cross_tab.n_local()); - if cross_tab.n_cols() <= cross_tab.n_rows() { - return (cross_tab, diagonal, indices); - } - diagonal.rotate_left(cross_tab.n_rows()); - indices.rotate_left(cross_tab.n_rows()); - let (c, ct) = cross_tab.into_parts(); - (CrossTab::with_transpose(ct, c), diagonal, indices) -} - -/// Bipartite SDDM in eliminated-major form; arrays are flat, split at `n_eliminated`. -#[derive(Clone)] -pub(crate) struct SddmMatrix { - pub(crate) cross_tab: CrossTab, - pub(crate) diagonal: Vec, - pub(crate) ground_edges: Vec, - pub(crate) grounding: Grounding, -} - -impl SddmMatrix { - pub(crate) fn n_eliminated(&self) -> usize { - self.cross_tab.n_rows() - } - - pub(crate) fn n_kept(&self) -> usize { - self.cross_tab.n_cols() - } - - pub(crate) fn diagonal_kept(&self) -> &[f64] { - &self.diagonal[self.n_eliminated()..] - } - - pub(crate) fn surplus_eliminated(&self) -> &[f64] { - &self.ground_edges[..self.n_eliminated()] - } - - pub(crate) fn surplus_kept(&self) -> &[f64] { - &self.ground_edges[self.n_eliminated()..] - } - - /// A floating matrix anchors one node instead, so it drops its surplus. - pub(super) fn float(&mut self) { - self.diagonal = magnitude_sums(&self.cross_tab); - self.ground_edges.fill(0.0); - self.grounding = Grounding::Floating; - } -} - -/// A connected component of a channel pair's cross-tab, converted to SDDM form. -#[derive(Clone)] -pub(crate) struct LocalComponent { - pub(crate) matrix: SddmMatrix, - pub(crate) form: MatrixForm, - pub(crate) coordinates: CoordinateMap, -} - -/// No diagonal scaling reaches weak dominance (Boman), so no SDDM form exists. -#[derive(Clone, Copy, Debug, PartialEq, Eq)] -pub(super) struct NotScalable; - -/// Evidence a dominance scaling exceeded tolerance or budget under [`ScalingFailure::Warn`]. -#[derive(Clone, Copy, Debug)] -pub(crate) struct UncertifiedScaling { - pub(crate) iterations: usize, - pub(crate) violation: f64, -} - -pub(super) fn convert( - cross_tab: CrossTab, - diagonal: Vec, - class: ComponentClass, - scaling: &ScalingConfig, -) -> Result<(LocalComponent, Option), NotScalable> { - match class { - ComponentClass::KnownLaplacian => { - convert_known_laplacian(cross_tab, diagonal).map(|component| (component, None)) - } - ComponentClass::General => convert_general(cross_tab, diagonal, scaling), - } -} - -/// Skipping the signed machinery keeps the dominant plain path at two streaming passes. -fn convert_known_laplacian( - cross_tab: CrossTab, - diagonal: Vec, -) -> Result { - let row_sums = adjacency_sums(&cross_tab); - let mut total_diagonal = 0.0; - let mut total_mismatch = 0.0; - for (&entry, &row_sum) in diagonal.iter().zip(row_sums.iter()) { - total_diagonal += entry; - total_mismatch += (entry - row_sum).abs(); - } - let n = cross_tab.n_local().max(1) as f64; - if total_mismatch > LAPLACIAN_VALIDATION_ULPS * f64::EPSILON * n.sqrt() * total_diagonal { - return Err(NotScalable); - } - Ok(LocalComponent { - matrix: SddmMatrix { - ground_edges: vec![0.0; cross_tab.n_local()], - cross_tab, - diagonal: row_sums, - grounding: Grounding::Floating, - }, - form: MatrixForm::Laplacian, - coordinates: CoordinateMap::default(), - }) -} - -fn convert_general( - cross_tab: CrossTab, - diagonal: Vec, - scaling: &ScalingConfig, -) -> Result<(LocalComponent, Option), NotScalable> { - let signs = folding_signs(&cross_tab); - let relaxation = dominance_scaling(&cross_tab, &diagonal, scaling)?; - if relaxation.violation > scaling.tolerance && scaling.on_failure == ScalingFailure::Error { - return Err(NotScalable); - } - let (component, clamped_deficit) = match signs { - Some(signs) => { - let factors: Vec = signs - .iter() - .zip(relaxation.scales.iter()) - .map(|(&sign, &scale)| sign * scale) - .collect(); - assemble(cross_tab, diagonal, factors, MatrixForm::Laplacian, scaling)? - } - None => { - let n_rows = cross_tab.n_rows(); - let mut factors = relaxation.scales; - for factor in &mut factors[n_rows..] { - *factor = -*factor; - } - assemble( - cross_tab, - diagonal, - factors, - MatrixForm::SignedPendingCover, - scaling, - )? - } - }; - let violation = relaxation.violation.max(clamped_deficit); - let uncertified = (violation > scaling.tolerance).then_some(UncertifiedScaling { - iterations: relaxation.iterations, - violation, - }); - Ok((component, uncertified)) -} - -/// A [`MatrixForm::Laplacian`] must fold to a Z-matrix, so a positive off-diagonal is an error. -fn assemble( - cross_tab: CrossTab, - diagonal: Vec, - factors: Vec, - form: MatrixForm, - scaling: &ScalingConfig, -) -> Result<(LocalComponent, f64), NotScalable> { - let n_rows = cross_tab.n_rows(); - let enforce_z = form == MatrixForm::Laplacian; - let (mut c, mut ct) = cross_tab.into_parts(); - let mut row_sums = Vec::with_capacity(c.nrows + ct.nrows); - for i in 0..n_rows { - let start = c.indptr[i] as usize; - let end = c.indptr[i + 1] as usize; - let columns = &c.indices[start..end]; - for (&j, value) in columns.iter().zip(&mut c.data[start..end]) { - let folded = -factors[i] * factors[n_rows + j as usize] * *value; - if !folded.is_finite() || (enforce_z && folded < 0.0) { - return Err(NotScalable); - } - *value = folded; - } - row_sums.push(c.data[start..end].iter().map(|v| v.abs()).sum()); - } - // Same product order as `c`, so each `Cᵀ` entry folds to the bits of its `C` twin. - for j in 0..ct.nrows { - let start = ct.indptr[j] as usize; - let end = ct.indptr[j + 1] as usize; - let rows = &ct.indices[start..end]; - for (&i, value) in rows.iter().zip(&mut ct.data[start..end]) { - *value *= -factors[i as usize] * factors[n_rows + j]; - } - row_sums.push(ct.data[start..end].iter().map(|v| v.abs()).sum()); - } - let cross_tab = CrossTab::with_transpose(c, ct); - - let scaled_diagonal: Vec = diagonal - .iter() - .zip(factors.iter()) - .map(|(&entry, &factor)| entry * factor * factor) - .collect(); - - // Canonical bipartite factors need no stored map: `fold` applies that sign flip by default. - let canonical = factors - .iter() - .enumerate() - .all(|(i, &factor)| factor == if i < n_rows { 1.0 } else { -1.0 }); - let coordinates = if canonical { - CoordinateMap::Canonical - } else { - CoordinateMap::Scaled(factors.into_boxed_slice()) - }; - - finalize( - cross_tab, - scaled_diagonal, - &row_sums, - coordinates, - form, - scaling, - ) -} - -/// Clamp roundoff deficits and retain surplus as ground edges. -fn finalize( - cross_tab: CrossTab, - mut scaled_diagonal: Vec, - row_sums: &[f64], - coordinates: CoordinateMap, - form: MatrixForm, - scaling: &ScalingConfig, -) -> Result<(LocalComponent, f64), NotScalable> { - let mut ground_edges = vec![0.0; cross_tab.n_local()]; - - let mut clamped_deficit = 0.0f64; - for ((diagonal, &row_sum), surplus) in scaled_diagonal - .iter_mut() - .zip(row_sums.iter()) - .zip(ground_edges.iter_mut()) - { - if !diagonal.is_finite() || *diagonal <= 0.0 { - return Err(NotScalable); - } - let deficit = (row_sum - *diagonal) / *diagonal; - if deficit > scaling.tolerance { - if scaling.on_failure == ScalingFailure::Error { - return Err(NotScalable); - } - clamped_deficit = clamped_deficit.max(deficit); - } - *diagonal = diagonal.max(row_sum); - *surplus = (*diagonal - row_sum).max(0.0); - } - Ok(( - LocalComponent { - matrix: SddmMatrix { - cross_tab, - diagonal: scaled_diagonal, - ground_edges, - grounding: Grounding::Grounded, - }, - form, - coordinates, - }, - clamped_deficit, - )) -} - -/// `None` when the component is frustrated: a negative cycle admits no signature. -fn folding_signs(cross_tab: &CrossTab) -> Option> { - let n = cross_tab.n_local(); - let mut signs = vec![0.0; n]; - if n == 0 { - return Some(signs); - } - for root in 0..n { - if signs[root] != 0.0 { - continue; - } - signs[root] = 1.0; - let mut stack = vec![root]; - while let Some(i) = stack.pop() { - for (j, value) in cross_tab.neighbors(i) { - let expected = -signs[i] * value.signum(); - if signs[j] == 0.0 { - signs[j] = expected; - stack.push(j); - } else if signs[j] != expected { - return None; - } - } - } - } - Some(signs) -} - -fn adjacency_sums(cross_tab: &CrossTab) -> Vec { - (0..cross_tab.n_local()) - .map(|i| cross_tab.neighbors(i).map(|(_, v)| v).sum()) - .collect() -} - -fn magnitude_sums(cross_tab: &CrossTab) -> Vec { - (0..cross_tab.n_local()) - .map(|i| cross_tab.neighbors(i).map(|(_, v)| v.abs()).sum()) - .collect() -} - -#[cfg(test)] -mod tests; diff --git a/crates/within/src/domain/factor_pairs/sddm/scaling.rs b/crates/within/src/domain/factor_pairs/sddm/scaling.rs deleted file mode 100644 index 7e8995dc..00000000 --- a/crates/within/src/domain/factor_pairs/sddm/scaling.rs +++ /dev/null @@ -1,401 +0,0 @@ -use super::NotScalable; -use crate::config::ScalingConfig; -use crate::csr_block::CsrBlock; -use crate::domain::CrossTab; -use crate::linalg::dot; - -pub(super) struct DominanceScaling { - pub(super) scales: Vec, - pub(super) iterations: usize, - /// Largest relative dominance violation at hand-over, clamped at 0. - pub(super) violation: f64, -} - -#[derive(Clone, Copy)] -enum BipartiteSide { - Rows, - Columns, -} - -struct NormalizedCrossOperator<'a> { - small_to_large: &'a CsrBlock, - large_to_small: &'a CsrBlock, - small_inv_sqrt: &'a [f64], - large_inv_sqrt: &'a [f64], - small_side: BipartiteSide, -} - -#[derive(Clone, Copy)] -enum CandidateKind { - Strict, - Boundary, -} - -impl<'a> NormalizedCrossOperator<'a> { - fn new(cross_tab: &'a CrossTab, inv_sqrt: &'a [f64]) -> Self { - let (row_inv_sqrt, col_inv_sqrt) = inv_sqrt.split_at(cross_tab.n_rows()); - let (small_to_large, large_to_small, small_inv_sqrt, large_inv_sqrt, small_side) = - if cross_tab.n_rows() <= cross_tab.n_cols() { - ( - cross_tab.ct(), - &cross_tab.c, - row_inv_sqrt, - col_inv_sqrt, - BipartiteSide::Rows, - ) - } else { - ( - &cross_tab.c, - cross_tab.ct(), - col_inv_sqrt, - row_inv_sqrt, - BipartiteSide::Columns, - ) - }; - Self { - small_to_large, - large_to_small, - small_inv_sqrt, - large_inv_sqrt, - small_side, - } - } - - fn small_len(&self) -> usize { - self.small_inv_sqrt.len() - } - - fn large_len(&self) -> usize { - self.large_inv_sqrt.len() - } - - fn apply_large(&self, input: &[f64], output: &mut [f64]) { - multiply_normalized( - self.small_to_large, - self.large_inv_sqrt, - self.small_inv_sqrt, - input, - output, - ); - } - - fn apply_small(&self, input: &[f64], output: &mut [f64]) { - multiply_normalized( - self.large_to_small, - self.small_inv_sqrt, - self.large_inv_sqrt, - input, - output, - ); - } - - fn apply_reduced(&self, input: &[f64], large_work: &mut [f64], output: &mut [f64]) { - // Eliminating the large side gives `I - MᵀM` on the cheaper side. - self.apply_large(input, large_work); - self.apply_small(large_work, output); - for (value, &input_value) in output.iter_mut().zip(input) { - *value = input_value - *value; - } - } - - fn reduced_rhs(&self, large_work: &mut [f64]) -> Vec { - large_work.fill(1.0); - let mut rhs = vec![0.0; self.small_len()]; - self.apply_small(large_work, &mut rhs); - for value in &mut rhs { - *value += 1.0; - } - rhs - } - - fn candidate_violation( - &self, - small: &[f64], - kind: CandidateKind, - large: &mut [f64], - small_image: &mut [f64], - ) -> Option { - if !small.iter().all(|value| value.is_finite() && *value > 0.0) { - return None; - } - self.apply_large(small, large); - if matches!(kind, CandidateKind::Strict) { - for value in large.iter_mut() { - *value += 1.0; - } - } - if !large.iter().all(|value| value.is_finite() && *value > 0.0) { - return None; - } - self.apply_small(large, small_image); - Some( - small_image - .iter() - .zip(small) - .map(|(&target, &value)| target / value - 1.0) - .fold(0.0f64, f64::max), - ) - } - - fn materialize_candidate( - &self, - small: &[f64], - large: &[f64], - violation: f64, - ) -> Result { - let mut mu: Vec = match self.small_side { - BipartiteSide::Rows => small.iter().chain(large).copied().collect(), - BipartiteSide::Columns => large.iter().chain(small).copied().collect(), - }; - let peak = mu.iter().copied().fold(0.0f64, f64::max); - if !peak.is_finite() || peak <= 0.0 { - return Err(NotScalable); - } - for value in mu.iter_mut() { - *value /= peak; - } - if !mu.iter().all(|value| value.is_finite() && *value > 0.0) { - return Err(NotScalable); - } - Ok(ScalingCandidate { mu, violation }) - } -} - -fn multiply_normalized( - matrix: &CsrBlock, - output_inv_sqrt: &[f64], - input_inv_sqrt: &[f64], - input: &[f64], - output: &mut [f64], -) { - debug_assert_eq!(matrix.nrows, output.len()); - debug_assert_eq!(matrix.ncols, input.len()); - for (i, value) in output.iter_mut().enumerate() { - *value = matrix - .row(i) - .map(|(j, entry)| entry.abs() * output_inv_sqrt[i] * input_inv_sqrt[j] * input[j]) - .sum(); - } -} - -struct ScalingCandidate { - mu: Vec, - violation: f64, -} - -struct ScalingResult { - candidate: ScalingCandidate, - iterations: usize, -} - -#[derive(Clone, Copy)] -enum ReducedCgStep { - Advanced, - ResidualExhausted, - NonPositiveCurvature, -} - -struct ReducedCg { - solution: Vec, - residual: Vec, - direction: Vec, - reduced_image: Vec, - residual_norm_squared: f64, -} - -impl ReducedCg { - fn new(rhs: Vec) -> Result { - let dimension = rhs.len(); - let residual_norm_squared = dot(&rhs, &rhs); - if !residual_norm_squared.is_finite() { - return Err(NotScalable); - } - Ok(Self { - solution: vec![0.0; dimension], - residual: rhs.clone(), - direction: rhs, - reduced_image: vec![0.0; dimension], - residual_norm_squared, - }) - } - - fn step( - &mut self, - operator: &NormalizedCrossOperator<'_>, - large_work: &mut [f64], - ) -> Result { - operator.apply_reduced(&self.direction, large_work, &mut self.reduced_image); - let direction_norm_squared = dot(&self.direction, &self.direction); - let curvature = dot(&self.direction, &self.reduced_image); - if !direction_norm_squared.is_finite() || !curvature.is_finite() { - return Err(NotScalable); - } - let slack = 64.0 * f64::EPSILON * direction_norm_squared; - if curvature <= slack { - return Ok(ReducedCgStep::NonPositiveCurvature); - } - - let alpha = self.residual_norm_squared / curvature; - if !alpha.is_finite() { - return Err(NotScalable); - } - for ((value, residual), (&direction, &image)) in self - .solution - .iter_mut() - .zip(&mut self.residual) - .zip(self.direction.iter().zip(&self.reduced_image)) - { - *value += alpha * direction; - *residual -= alpha * image; - } - - let next_residual_norm_squared = dot(&self.residual, &self.residual); - if !next_residual_norm_squared.is_finite() { - return Err(NotScalable); - } - if next_residual_norm_squared == 0.0 { - return Ok(ReducedCgStep::ResidualExhausted); - } - let beta = next_residual_norm_squared / self.residual_norm_squared; - for (direction, &residual) in self.direction.iter_mut().zip(&self.residual) { - *direction = residual + beta * *direction; - } - self.residual_norm_squared = next_residual_norm_squared; - Ok(ReducedCgStep::Advanced) - } -} - -fn reduced_cg_scaling( - operator: &NormalizedCrossOperator<'_>, - initial: ScalingCandidate, - scaling: &ScalingConfig, -) -> Result { - let mut large_work = vec![0.0; operator.large_len()]; - let rhs = operator.reduced_rhs(&mut large_work); - let mut cg = ReducedCg::new(rhs)?; - let mut best = initial; - let mut small_image = vec![0.0; operator.small_len()]; - - for iteration in 1..=scaling.max_iterations { - let step = cg.step(operator, &mut large_work)?; - if matches!(step, ReducedCgStep::NonPositiveCurvature) { - // A singular boundary exposes its Perron vector as the zero-curvature direction. - let boundary_small: Vec = cg.direction.iter().map(|value| value.abs()).collect(); - if let Some(violation) = operator.candidate_violation( - &boundary_small, - CandidateKind::Boundary, - &mut large_work, - &mut small_image, - ) { - if violation < best.violation { - if let Ok(candidate) = - operator.materialize_candidate(&boundary_small, &large_work, violation) - { - best = candidate; - } - } - } - return Ok(ScalingResult { - candidate: best, - iterations: iteration, - }); - } - - // A certificate check costs two extra passes; past the cheap early exits, probe on a stride. - let last = - matches!(step, ReducedCgStep::ResidualExhausted) || iteration == scaling.max_iterations; - if last || iteration <= 8 || iteration % 4 == 0 { - if let Some(violation) = operator.candidate_violation( - &cg.solution, - CandidateKind::Strict, - &mut large_work, - &mut small_image, - ) { - if violation < best.violation { - if let Ok(candidate) = - operator.materialize_candidate(&cg.solution, &large_work, violation) - { - best = candidate; - } - } - if best.violation <= scaling.tolerance { - return Ok(ScalingResult { - candidate: best, - iterations: iteration, - }); - } - } - } - - if matches!(step, ReducedCgStep::ResidualExhausted) { - return Ok(ScalingResult { - candidate: best, - iterations: iteration, - }); - } - } - Ok(ScalingResult { - candidate: best, - iterations: scaling.max_iterations, - }) -} - -pub(super) fn dominance_scaling( - cross_tab: &CrossTab, - diagonal: &[f64], - scaling: &ScalingConfig, -) -> Result { - let n = cross_tab.n_local(); - debug_assert_eq!(diagonal.len(), n); - if diagonal.iter().any(|d| !d.is_finite() || *d <= 0.0) { - return Err(NotScalable); - } - - let already_dominant = (0..n).all(|i| { - let row_sum: f64 = cross_tab.neighbors(i).map(|(_, value)| value.abs()).sum(); - row_sum <= diagonal[i] * (1.0 + scaling.tolerance) - }); - if already_dominant { - return Ok(DominanceScaling { - scales: vec![1.0; n], - iterations: 0, - violation: 0.0, - }); - } - - let inv_sqrt: Vec = diagonal.iter().map(|d| 1.0 / d.sqrt()).collect(); - let violation_of = |mu: &[f64]| { - (0..n) - .map(|i| { - let target: f64 = cross_tab - .neighbors(i) - .map(|(j, value)| value.abs() * inv_sqrt[i] * inv_sqrt[j] * mu[j]) - .sum(); - target / mu[i] - 1.0 - }) - .fold(0.0f64, f64::max) - }; - - let initial_mu = vec![1.0; n]; - let initial = ScalingCandidate { - violation: violation_of(&initial_mu), - mu: initial_mu, - }; - let operator = NormalizedCrossOperator::new(cross_tab, &inv_sqrt); - let result = reduced_cg_scaling(&operator, initial, scaling)?; - let scales: Vec = result - .candidate - .mu - .iter() - .zip(inv_sqrt.iter()) - .map(|(&value, &normalizer)| value * normalizer) - .collect(); - if !scales.iter().all(|value| value.is_finite() && *value > 0.0) { - return Err(NotScalable); - } - Ok(DominanceScaling { - scales, - iterations: result.iterations, - violation: result.candidate.violation.max(0.0), - }) -} diff --git a/crates/within/src/domain/factor_pairs/sddm/tests.rs b/crates/within/src/domain/factor_pairs/sddm/tests.rs deleted file mode 100644 index d1905dd3..00000000 --- a/crates/within/src/domain/factor_pairs/sddm/tests.rs +++ /dev/null @@ -1,340 +0,0 @@ -use super::*; -use crate::csr_block::CsrBlock; - -impl LocalComponent { - pub(crate) fn plain_for_test(cross_tab: CrossTab, diagonal: Vec) -> Self { - convert( - cross_tab, - diagonal, - ComponentClass::KnownLaplacian, - &ScalingConfig::default(), - ) - .expect("plain test component must convert to SDDM") - .0 - } - - pub(crate) fn general_for_test(cross_tab: CrossTab, diagonal: Vec) -> Self { - let globals = (0..cross_tab.n_local() as u32).collect(); - let (cross_tab, diagonal, _) = super::orient_for_elimination(cross_tab, diagonal, globals); - convert( - cross_tab, - diagonal, - ComponentClass::General, - &ScalingConfig::default(), - ) - .expect("general test component must convert to SDDM") - .0 - } - - /// Assemble under externally supplied congruence factors: the relaxation's - /// certified scaling is valid but not unique, so tests that construct a - /// component backward from known factors pin those factors here. - pub(crate) fn with_factors_for_test( - cross_tab: CrossTab, - diagonal: Vec, - factors: &[f64], - ) -> Self { - assemble( - cross_tab, - diagonal, - factors.to_vec(), - MatrixForm::Laplacian, - &ScalingConfig::default(), - ) - .expect("test factors must fold to SDDM") - .0 - } -} - -impl SddmMatrix { - pub(crate) fn from_dense_for_test( - table: &[f64], - n_rows: usize, - n_cols: usize, - diagonal: Vec, - grounding: Grounding, - ) -> Self { - let cross_tab = CrossTab::from_dense_for_test(table, n_rows, n_cols); - debug_assert_eq!(diagonal.len(), cross_tab.n_local()); - let ground_edges = adjacency_sums(&cross_tab) - .iter() - .zip(&diagonal) - .map(|(sum, d)| (d - sum).max(0.0)) - .collect(); - Self { - cross_tab, - diagonal, - ground_edges, - grounding, - } - } - - pub(crate) fn laplacian_for_test(table: &[f64], n_rows: usize, n_cols: usize) -> Self { - let diagonal = adjacency_sums(&CrossTab::from_dense_for_test(table, n_rows, n_cols)); - Self::from_dense_for_test(table, n_rows, n_cols, diagonal, Grounding::Floating) - } -} - -/// Chain: row k joins cols k and k+1; one mid-chain deficit makes certification need O(m²) iterations. -fn chain_with_mid_deficit(m: usize) -> (CrossTab, Vec) { - let c = CsrBlock { - indptr: (0..=m as u32).map(|k| 2 * k).collect(), - indices: (0..m as u32).flat_map(|k| [k, k + 1]).collect(), - data: vec![1.0; 2 * m], - nrows: m, - ncols: m + 1, - }; - let mut diagonal = vec![2.0; 2 * m + 1]; - diagonal[m + m / 2] = 2.0 * (1.0 - 1e-3); - (CrossTab::eager(c), diagonal) -} - -#[test] -fn long_feasible_chain_certifies_with_reduced_cg() { - let (cross_tab, diagonal) = chain_with_mid_deficit(200); - let scaling = ScalingConfig { - max_iterations: 128, - ..Default::default() - }; - let result = dominance_scaling(&cross_tab, &diagonal, &scaling).expect("chain scales"); - assert!( - result.violation <= scaling.tolerance, - "chain must certify, got violation {:.3e}", - result.violation, - ); - assert!( - result.iterations < scaling.max_iterations, - "reduced CG exhausted the {} iteration budget", - scaling.max_iterations, - ); -} - -#[test] -fn accelerating_convergence_case_certifies_without_extrapolation() { - let table = [ - 0.04023094100803415, - 0.9632989615043355, - 0.003631923482738449, - 0.0015750938044434863, - 0.0006618152745868568, - 0.039346817354759124, - ]; - let cross_tab = CrossTab::from_dense_for_test(&table, 3, 2); - let scaling = ScalingConfig { - max_iterations: 16, - ..Default::default() - }; - let result = dominance_scaling(&cross_tab, &[1.0; 5], &scaling).expect("matrix scales"); - assert!( - result.violation <= scaling.tolerance, - "matrix must certify, got violation {:.3e} after {} iterations", - result.violation, - result.iterations, - ); - assert!(result.iterations <= 2); -} - -#[test] -fn indefinite_comparison_matrix_remains_uncertified() { - let cross_tab = CrossTab::from_dense_for_test(&[1.0, -1.0, 2.0, -2.0], 2, 2); - let diagonal: [f64; 4] = [1.0, 2.0, 1.0, 2.0]; - let scaling = ScalingConfig::default(); - let result = dominance_scaling(&cross_tab, &diagonal, &scaling).unwrap(); - assert!(result.violation > scaling.tolerance); - assert!(result.iterations < scaling.max_iterations); -} - -fn assert_sddm(component: &LocalComponent) { - assert!(component - .matrix - .cross_tab - .c - .data - .iter() - .all(|value| *value >= 0.0)); - let sums = adjacency_sums(&component.matrix.cross_tab); - for ((&diagonal, &row_sum), &surplus) in component - .matrix - .diagonal - .iter() - .zip(sums.iter()) - .zip(component.matrix.ground_edges.iter()) - { - assert!(diagonal >= row_sum); - assert!((diagonal - row_sum - surplus).abs() <= 1e-12 * diagonal); - } -} - -#[test] -fn known_laplacian_has_canonical_coordinates_and_no_ground() { - let (component, uncertified) = convert( - CrossTab::from_dense_for_test(&[2.0, 1.0, 0.0, 3.0], 2, 2), - vec![3.0, 3.0, 2.0, 4.0], - ComponentClass::KnownLaplacian, - &ScalingConfig::default(), - ) - .unwrap(); - assert_eq!(component.form, MatrixForm::Laplacian); - assert_eq!(component.matrix.grounding, Grounding::Floating); - assert!(matches!(component.coordinates, CoordinateMap::Canonical)); - assert!(uncertified.is_none()); - assert_sddm(&component); -} - -#[test] -fn known_laplacian_claim_is_checked() { - // Structural surplus contradicts the Laplacian claim, so it must fail loudly. - let result = convert( - CrossTab::from_dense_for_test(&[2.0, 1.0, 0.0, 3.0], 2, 2), - vec![4.0, 3.0, 2.0, 4.0], - ComponentClass::KnownLaplacian, - &ScalingConfig::default(), - ); - assert!(matches!(result, Err(NotScalable))); -} - -#[test] -fn frustrated_component_stores_single_signed_operator() { - let (component, uncertified) = convert( - CrossTab::from_dense_for_test(&[1.0, 1.0, 1.0, -1.0], 2, 2), - vec![2.0, 2.0, 2.0, 2.0], - ComponentClass::General, - &ScalingConfig::default(), - ) - .unwrap(); - assert!(uncertified.is_none()); - // The matrix stays single-sized and signed; the cover is deferred to factor time. - assert_eq!(component.form, MatrixForm::SignedPendingCover); - assert_eq!(component.matrix.cross_tab.c.nrows, 2); - assert_eq!(component.matrix.cross_tab.c.ncols, 2); - let mut dense = [[0.0; 2]; 2]; - for (i, row) in dense.iter_mut().enumerate() { - for (j, value) in component.matrix.cross_tab.c.row(i) { - row[j] = value; - } - } - assert_eq!(dense, [[1.0, 1.0], [1.0, -1.0]]); - // Zero surplus and the canonical bipartite sign flip, so no explicit factor map is stored. - assert!(component.matrix.ground_edges.iter().all(|&s| s == 0.0)); - assert!(matches!(component.coordinates, CoordinateMap::Canonical)); -} - -#[test] -fn scalable_signed_component_produces_valid_grounded_sddm() { - let d = [2.0, -0.5, -1.0, 4.0, 0.25]; - let c_hat = [[1.0, 2.0, 0.5], [3.0, 0.0, 1.5]]; - let diag_hat = [4.0, 5.0, 4.5, 2.5, 2.375]; - let mut raw = vec![0.0; 6]; - for i in 0..2 { - for j in 0..3 { - raw[i * 3 + j] = c_hat[i][j] / (d[i] * d[2 + j]); - } - } - let (component, uncertified) = convert( - CrossTab::from_dense_for_test(&raw, 2, 3), - (0..5).map(|i| diag_hat[i] / (d[i] * d[i])).collect(), - ComponentClass::General, - &ScalingConfig::default(), - ) - .unwrap(); - assert_eq!(component.form, MatrixForm::Laplacian); - assert_eq!(component.matrix.grounding, Grounding::Grounded); - assert!(uncertified.is_none()); - assert_sddm(&component); -} - -#[test] -fn non_scalable_component_errors_under_error_mode() { - let result = convert( - CrossTab::from_dense_for_test(&[1.0, -1.0, 2.0, -2.0], 2, 2), - vec![1.0, 2.0, 1.0, 2.0], - ComponentClass::General, - &ScalingConfig { - on_failure: ScalingFailure::Error, - ..Default::default() - }, - ); - assert!(matches!(result, Err(NotScalable))); -} - -#[test] -fn non_scalable_component_warns_and_clamps_under_warn_mode() { - let config = ScalingConfig { - max_iterations: 32, - ..Default::default() - }; - assert_eq!(config.on_failure, ScalingFailure::Warn); - let (component, uncertified) = convert( - CrossTab::from_dense_for_test(&[1.0, -1.0, 2.0, -2.0], 2, 2), - vec![1.0, 2.0, 1.0, 2.0], - ComponentClass::General, - &config, - ) - .unwrap(); - let uncertified = uncertified.expect("no dominant scaling exists"); - assert!(uncertified.violation > config.tolerance); - // Clamping restores the SDDM invariant, so the hand-over is usable. - assert_sddm(&component); -} - -#[test] -fn barely_pd_surplus_is_structural() { - let surplus = 5e-10; - let (component, _) = convert( - CrossTab::from_dense_for_test(&[1.0], 1, 1), - vec![1.0 + surplus, 1.0], - ComponentClass::General, - &ScalingConfig::default(), - ) - .unwrap(); - assert_eq!(component.form, MatrixForm::Laplacian); - assert_eq!(component.matrix.grounding, Grounding::Grounded); - assert!((component.matrix.ground_edges[0] - surplus).abs() < 1e-15); - assert_sddm(&component); -} - -#[test] -fn long_weighted_runs_preserve_known_laplacian_balance() { - use crate::channel::{Channel, ChannelPair}; - use crate::domain::{Design, Effect, PreparedDesign}; - let n = 1_000_000; - // The first case exposed diagonal drift; the second also has long cell sums. - for [na, nb, nc] in [[317, 101, 7], [5, 2, 3]] { - let a: Vec = (0..n).map(|i| (i % na) as u32).collect(); - let b: Vec = (0..n).map(|i| ((i / 3 + i % 7) % nb) as u32).collect(); - let c: Vec = (0..n).map(|i| ((i / 19) % nc) as u32).collect(); - let weights: Vec = (0..n).map(|i| 0.7 + (i % 13) as f64 / 7.).collect(); - let effects = [&a, &b, &c].map(|v| Effect::new(v, true, []).unwrap()); - let prepared = PreparedDesign::new(Design::new(effects).unwrap(), Some(&weights)).unwrap(); - let pair = ChannelPair { - rows: Channel { term: 1, column: 0 }, - cols: Channel { term: 2, column: 0 }, - }; - let (cross_tab, _) = CrossTab::build_for_pair(&prepared, pair); - let diagonal = [ - prepared.channel_diagonal(pair.rows), - prepared.channel_diagonal(pair.cols), - ] - .concat(); - assert!(convert_known_laplacian(cross_tab, diagonal).is_ok()); - // Integer multiplicities provide an independent short-sum oracle; no - // long observation fold or production compensation is reused here. - for (channel, levels, width) in [(pair.rows, &b, nb), (pair.cols, &c, nc)] { - let mut counts = vec![[0u32; 13]; width]; - for (i, &level) in levels.iter().enumerate() { - counts[level as usize][i % 13] += 1; - } - for (&actual, count) in prepared.channel_diagonal(channel).iter().zip(counts) { - let expected: f64 = count - .into_iter() - .enumerate() - .map(|(r, count)| { - let s = (0.7 + r as f64 / 7.).sqrt(); - f64::from(count) * (s * s) - }) - .sum(); - assert!((actual - expected).abs() <= 16.0 * f64::EPSILON * expected); - } - } - } -} diff --git a/crates/within/src/domain/prepared.rs b/crates/within/src/domain/prepared.rs index 6646399e..45fd60f8 100644 --- a/crates/within/src/domain/prepared.rs +++ b/crates/within/src/domain/prepared.rs @@ -3,7 +3,7 @@ use std::sync::OnceLock; use rayon::prelude::*; use super::{row_weight, Column, Design, ScatterSplit, Term, TermReparam}; -use crate::channel::{Channel, CoefficientPosition}; +use crate::channel::CoefficientPosition; use crate::linalg::add_compensated; use crate::BuildError; @@ -139,13 +139,6 @@ impl<'a> PreparedDesign<'a> { }) } - /// Channel `channel`'s `n_levels` entries of `diag(AᵀA)`. - pub(crate) fn channel_diagonal(&self, channel: Channel) -> &[f64] { - let term = &self.design.terms[channel.term]; - let dofs = term.column_dofs(channel.column); - &self.diagonal(channel.term)[dofs.start - term.offset..dofs.end - term.offset] - } - /// `diag(AᵀA)` in solve coordinates. pub(crate) fn gram_diagonal(&self) -> Vec { let blocks: Vec<&[f64]> = (0..self.design.n_factors()) diff --git a/crates/within/src/error.rs b/crates/within/src/error.rs index 09b74453..0ebd4f46 100644 --- a/crates/within/src/error.rs +++ b/crates/within/src/error.rs @@ -96,8 +96,8 @@ pub enum BuildError { /// Actual column count of the supplied preconditioner. actual_cols: usize, }, - /// Every dominance comparison is `>`, so a NaN or infinite slack silently certifies anything. - #[error("scaling tolerance must be finite and non-negative, got {value}")] + /// A scaling tolerance outside `[0, 1)`. + #[error("scaling tolerance must lie in [0, 1), got {value}")] InvalidScalingTolerance { /// The offending value. value: f64, diff --git a/crates/within/src/lib.rs b/crates/within/src/lib.rs index cfe41c56..90a432fb 100644 --- a/crates/within/src/lib.rs +++ b/crates/within/src/lib.rs @@ -26,9 +26,8 @@ pub mod config; pub mod error; -pub(crate) mod block_elim; +pub(crate) mod block_solver; pub(crate) mod channel; -pub(crate) mod csr_block; pub(crate) mod domain; pub(crate) mod linalg; pub(crate) mod operator; diff --git a/crates/within/src/linalg.rs b/crates/within/src/linalg.rs index 9566d4e3..a4219789 100644 --- a/crates/within/src/linalg.rs +++ b/crates/within/src/linalg.rs @@ -13,14 +13,6 @@ pub(crate) fn add_compensated(sum: &mut f64, correction: &mut f64, value: f64) { *sum = next; } -pub(crate) fn compensated_sum(values: &[f64]) -> f64 { - let (mut sum, mut correction) = (0.0, 0.0); - for &value in values { - add_compensated(&mut sum, &mut correction, value); - } - sum + correction -} - pub(crate) fn dot(left: &[f64], right: &[f64]) -> f64 { left.iter().zip(right).map(|(&x, &y)| x * y).sum() } diff --git a/crates/within/src/operator/schwarz.rs b/crates/within/src/operator/schwarz.rs index 39cfd393..ecd1f165 100644 --- a/crates/within/src/operator/schwarz.rs +++ b/crates/within/src/operator/schwarz.rs @@ -7,7 +7,7 @@ use serde::{Deserialize, Serialize}; use std::sync::Arc; use std::time::{Duration, Instant}; -use crate::block_elim::BlockElimSolver; +use crate::block_solver::BlockSolver; use crate::config::{LocalSolverConfig, PreconditionerConfig, Staleness}; use crate::domain::{LocalDomain, PreparedDesign}; use crate::{BuildError, BuildWarning}; @@ -18,7 +18,7 @@ mod tests; /// Concrete additive Schwarz type used in the parent crate. #[derive(Clone, Serialize, Deserialize)] pub(crate) struct FeSchwarz { - inner: SchwarzPreconditioner, + inner: SchwarzPreconditioner, local_solver: LocalSolverConfig, } @@ -97,7 +97,7 @@ pub(crate) fn build_additive( ) -> Result { let entries = domains .into_par_iter() - .map(|domain| build_entry(domain, local_solver)) + .map(build_entry) .collect::, BuildError>>()?; Ok(FeSchwarz { inner: SchwarzPreconditioner::with_n_dofs(entries, n_dofs), @@ -105,13 +105,8 @@ pub(crate) fn build_additive( }) } -/// Build a single `SubdomainEntry` from a pre-built CrossTab. -pub(crate) fn build_entry( - domain: LocalDomain, - config: &LocalSolverConfig, -) -> Result, BuildError> { - let LocalDomain { core, component } = domain; - let solver = BlockElimSolver::build(component, config)?; +pub(crate) fn build_entry(domain: LocalDomain) -> Result, BuildError> { + let LocalDomain { core, solver } = domain; SubdomainEntry::try_new(core, solver).map_err(BuildError::Preconditioner) } diff --git a/crates/within/src/operator/schwarz/tests.rs b/crates/within/src/operator/schwarz/tests.rs index 6c632d80..c6f74a4f 100644 --- a/crates/within/src/operator/schwarz/tests.rs +++ b/crates/within/src/operator/schwarz/tests.rs @@ -7,11 +7,7 @@ use crate::config::{ ApproxCholConfig, ApproxSchurConfig, LocalSolverConfig, SchurMode, DEFAULT_DENSE_SCHUR_THRESHOLD, }; -use schwarz_precond::SubdomainCore; - -use crate::csr_block::CsrBlock; -use crate::domain::{build_local_domains, LocalDomain, PreparedDesign}; -use crate::domain::{CrossTab, LocalComponent}; +use crate::domain::{build_local_domains, Design, LocalDomain, PreparedDesign}; use crate::operator::schwarz::build_additive; use schwarz_precond::{LocalSolver, Operator}; @@ -25,92 +21,42 @@ fn make_test_data() -> (PreparedDesign<'static>, Vec) { (design, domain_pairs) } -fn synthetic_sparse_cross_tab(n_keep: usize, elim_ratio: usize) -> (CrossTab, Vec) { - let n_rows = n_keep * elim_ratio; - let n_cols = n_keep; - let mut indptr = Vec::with_capacity(n_rows + 1); - let mut indices = Vec::with_capacity(n_rows * 3); - let mut data = Vec::with_capacity(n_rows * 3); - let mut row_diag = vec![0.0; n_rows]; - let mut col_diag = vec![0.0; n_cols]; - - indptr.push(0); - for (i, row_diag_i) in row_diag.iter_mut().enumerate().take(n_rows) { - let mut row = [ - (i % n_cols, 1.0), - ((i + 1) % n_cols, 0.8), - ((i.wrapping_mul(17).wrapping_add(3)) % n_cols, 0.6), - ]; - row.sort_unstable_by_key(|&(col, _)| col); - - let mut row_sum = 0.0; - let mut cursor = 0usize; - while cursor < row.len() { - let col = row[cursor].0; - let mut value = row[cursor].1; - cursor += 1; - while cursor < row.len() && row[cursor].0 == col { - value += row[cursor].1; - cursor += 1; - } - - indices.push(col as u32); - data.push(value); - row_sum += value; - col_diag[col] += value; +/// One plain pair: `n_keep · elim_ratio` levels, each in three weighted cells of `n_keep`. +fn synthetic_star_design(n_keep: usize, elim_ratio: usize) -> PreparedDesign<'static> { + let (mut rows, mut cols, mut weights) = (Vec::new(), Vec::new(), Vec::new()); + for i in 0..n_keep * elim_ratio { + for (col, w) in [ + (i % n_keep, 1.0), + ((i + 1) % n_keep, 0.8), + ((i.wrapping_mul(17).wrapping_add(3)) % n_keep, 0.6), + ] { + rows.push(i as u32); + cols.push(col as u32); + weights.push(w); } - - *row_diag_i = row_sum; - indptr.push(indices.len() as u32); } - - let c = CsrBlock { - indptr, - indices, - data, - nrows: n_rows, - ncols: n_cols, - }; - ( - CrossTab::eager(c), - row_diag.into_iter().chain(col_diag).collect(), - ) -} - -fn make_nested_block_elim_domain_pairs( - n_keep: usize, - elim_ratio: usize, - n_subdomains: usize, -) -> (usize, Vec) { - let (cross_tab, block_diagonals) = synthetic_sparse_cross_tab(n_keep, elim_ratio); - let n_local = cross_tab.n_local(); - let global_indices: Vec = (0..n_local as u32).collect(); - - let domain_pairs = (0..n_subdomains) - .map(|_| LocalDomain { - core: SubdomainCore::uniform(global_indices.clone()), - component: LocalComponent::plain_for_test(cross_tab.clone(), block_diagonals.clone()), - }) - .collect(); - (n_local, domain_pairs) + let design = Design::from_levels_for_test(vec![rows, cols]); + PreparedDesign::new(design, Some(&weights)).expect("valid weights") } fn run_block_elim_nested_rayon_case() { let n_keep = 1_024; let elim_ratio = 32; - let (n_dofs, domain_pairs) = make_nested_block_elim_domain_pairs(n_keep, elim_ratio, 2); let config = LocalSolverConfig { approx_chol: ApproxCholConfig { split_merge: Some(8), seed: 42, }, - schur: SchurMode::Approximate(ApproxSchurConfig { - seed: 7, - ..Default::default() - }), dense_threshold: DEFAULT_DENSE_SCHUR_THRESHOLD, ..Default::default() }; + let design = synthetic_star_design(n_keep, elim_ratio); + let n_dofs = design.design.n_dofs; + let (mut domain_pairs, _) = build_local_domains(&design, &config).expect("domains"); + let [domain] = &domain_pairs[..] else { + panic!("one component expected, got {}", domain_pairs.len()) + }; + domain_pairs.push(domain.clone()); let rhs: Vec = (0..n_dofs).map(|i| ((i % 29) as f64) - 14.0).collect(); let schwarz = @@ -197,20 +143,25 @@ fn test_build_additive() { /// Three-entry eliminated stars, so clique sampling is not /// deterministic-exact and the two routes are separable. fn small_subdomain_solve(schur: SchurMode, dense_threshold: usize) -> Vec { - let (cross_tab, block_diagonals) = synthetic_sparse_cross_tab(8, 4); - let component = LocalComponent::plain_for_test(cross_tab, block_diagonals); let config = LocalSolverConfig { approx_chol: ApproxCholConfig::default(), schur, dense_threshold, ..Default::default() }; - let solver = - crate::block_elim::BlockElimSolver::build(component, &config).expect("block-elim build"); - let r: Vec = (0..solver.n_local()) + let (domains, _) = build_local_domains(&synthetic_star_design(8, 4), &config).expect("domains"); + let [domain] = &domains[..] else { + panic!("one component expected, got {}", domains.len()) + }; + let mut x: Vec = (0..domain.solver.n_local()) .map(|i| if i % 2 == 0 { 1.0 } else { -1.0 }) .collect(); - solver.solve_for_test(&r, 0.0) + let mut workspace = vec![0.0; domain.solver.workspace_len()]; + domain + .solver + .solve_local(&mut x, &mut workspace, false) + .expect("local solve"); + x } #[test] @@ -218,10 +169,7 @@ fn test_schur_mode_is_outranked_by_the_exact_route_for_small_reduced_system() { // Agreeing is only possible if both took the exact route. let exact = small_subdomain_solve(SchurMode::Exact, DEFAULT_DENSE_SCHUR_THRESHOLD); let sampled_mode = small_subdomain_solve( - SchurMode::Approximate(ApproxSchurConfig { - seed: 7, - ..Default::default() - }), + SchurMode::Approximate(ApproxSchurConfig::default()), DEFAULT_DENSE_SCHUR_THRESHOLD, ); assert_eq!(exact, sampled_mode); @@ -230,13 +178,8 @@ fn test_schur_mode_is_outranked_by_the_exact_route_for_small_reduced_system() { #[test] fn test_dense_threshold_zero_disables_exact_route() { let exact = small_subdomain_solve(SchurMode::Exact, DEFAULT_DENSE_SCHUR_THRESHOLD); - let approximate = small_subdomain_solve( - SchurMode::Approximate(ApproxSchurConfig { - seed: 7, - ..Default::default() - }), - 0, - ); + let approximate = + small_subdomain_solve(SchurMode::Approximate(ApproxSchurConfig::default()), 0); assert_ne!(exact, approximate); } diff --git a/crates/within/src/solver/tests.rs b/crates/within/src/solver/tests.rs index 98b4b50b..a99e2333 100644 --- a/crates/within/src/solver/tests.rs +++ b/crates/within/src/solver/tests.rs @@ -4,8 +4,8 @@ use rstest::rstest; use super::{CoefficientAddress, CoefficientLayout}; use crate::channel::Channel; -use crate::config::{LocalSolverConfig, LsmrOptions, DEFAULT_DENSE_SCHUR_THRESHOLD}; -use crate::domain::{build_local_domains, Design, Grounding, MatrixForm, PreparedDesign}; +use crate::config::{LocalSolverConfig, LsmrOptions}; +use crate::domain::Design; use crate::{Effect, PreconditionerConfig, Solver}; /// DGP kept in lockstep with `surplus_component_sampled_matches_exact_reduction` @@ -29,27 +29,6 @@ fn positive_slope_only_panel() -> (Vec, Vec, Vec) { (f, g, z) } -#[test] -fn positive_slope_only_pair_grounds_beyond_dense_threshold() { - let (f, g, z) = positive_slope_only_panel(); - let effects = vec![ - Effect::new(&f, false, [&z[..]]).expect("slope effect"), - Effect::new(&g, true, []).expect("plain effect"), - ]; - let prepared = PreparedDesign::unweighted_for_test(Design::new(effects).expect("design")); - let (domains, warnings) = - build_local_domains(&prepared, &LocalSolverConfig::default()).expect("domains"); - assert!( - domains.iter().any(|ld| { - let ct = &ld.component.matrix.cross_tab; - ld.component.form == MatrixForm::Laplacian - && ld.component.matrix.grounding == Grounding::Grounded - && ct.n_rows().min(ct.n_cols()) > DEFAULT_DENSE_SCHUR_THRESHOLD - }), - "fixture must ground a component past the dense threshold (warnings: {warnings:?})" - ); -} - #[test] fn an_alias_its_scaling_misses_keeps_its_surplus() { // The null vector spans 1e200, beyond the relaxation; floating would drop a third of the diagonal. diff --git a/crates/within/tests/error_paths.rs b/crates/within/tests/error_paths.rs index 8e5c721b..6998feb8 100644 --- a/crates/within/tests/error_paths.rs +++ b/crates/within/tests/error_paths.rs @@ -211,7 +211,7 @@ fn test_within_error_source_chains_through_transparent_wrapper() { /// silently — no `UnscalableComponent` under `Error`, and no `BuildWarning` under `Warn` either. #[rstest] fn test_invalid_local_solver_rejected( - #[values(-1e-9, f64::NAN, f64::INFINITY)] value: f64, + #[values(-1e-9, 1.0, f64::NAN, f64::INFINITY)] value: f64, #[values(false, true)] adaptive: bool, #[values(1, 2)] n_factors: usize, ) { diff --git a/crates/within/tests/fixtures/preconditioner_v21.postcard b/crates/within/tests/fixtures/preconditioner_v21.postcard index ee2f4bb5..be7e672f 100644 Binary files a/crates/within/tests/fixtures/preconditioner_v21.postcard and b/crates/within/tests/fixtures/preconditioner_v21.postcard differ diff --git a/crates/within/tests/solver.rs b/crates/within/tests/solver.rs index 24c0f4e5..59709dc3 100644 --- a/crates/within/tests/solver.rs +++ b/crates/within/tests/solver.rs @@ -193,7 +193,7 @@ fn fully_specified_additive() -> PreconditionerConfig { seed: 41, split_merge: Some(3), }, - schur: SchurMode::Approximate(ApproxSchurConfig { seed: 17, split: 2 }), + schur: SchurMode::Approximate(ApproxSchurConfig { split: 2 }), dense_threshold: 0, scaling: ScalingConfig { tolerance: 1e-7, diff --git a/crates/within/tests/warnings.rs b/crates/within/tests/warnings.rs index 615e6e59..aadcc287 100644 --- a/crates/within/tests/warnings.rs +++ b/crates/within/tests/warnings.rs @@ -6,14 +6,14 @@ use within::{ ScalingConfig, ScalingFailure, Solver, }; -// Two crossed slope-only factors: their signed slope cross-block is not -// diagonally dominant. A zero-tolerance / zero-iteration scaling policy reports +// Two crossed slope-only factors: level a0's cross entries sum to 4 > its diagonal 2, +// so the block needs a dominance scaling. A zero-tolerance / zero-iteration policy reports // that as an uncertified-scaling warning deterministically, without relying on // the (deliberately rare) default-tolerance path. const LA: [u32; 4] = [0, 0, 1, 1]; const LB: [u32; 4] = [0, 1, 0, 1]; const ZA: [f64; 4] = [1.0, 1.0, 1.0, 1.0]; -const ZB: [f64; 4] = [1.0, -1.0, 2.0, -2.0]; +const ZB: [f64; 4] = [3.0, -1.0, 2.0, -4.0]; fn warning_effects() -> Vec> { vec![ @@ -22,19 +22,6 @@ fn warning_effects() -> Vec> { ] } -fn strict_scaling() -> PreconditionerConfig { - PreconditionerConfig::Additive { - local_solver: LocalSolverConfig { - scaling: ScalingConfig { - tolerance: 0.0, - max_iterations: 0, - on_failure: ScalingFailure::Warn, - }, - ..LocalSolverConfig::default() - }, - } -} - fn has_scaling_warning(warnings: &[BuildWarning]) -> bool { warnings .iter() @@ -44,7 +31,7 @@ fn has_scaling_warning(warnings: &[BuildWarning]) -> bool { #[test] fn free_solve_surfaces_build_warnings() { let y = [1.0, 2.0, 0.5, -1.0]; - let cfg = strict_scaling(); + let cfg = scaling(0.0, ScalingFailure::Warn); let result = solve(warning_effects(), &y, None, None, &cfg).expect("solve"); assert!( @@ -63,7 +50,7 @@ fn free_solve_batch_surfaces_build_warnings() { let y0 = [1.0, 2.0, 0.5, -1.0]; let y1 = [0.2, -0.7, 1.3, 0.4]; let ys: [&[f64]; 2] = [&y0, &y1]; - let cfg = strict_scaling(); + let cfg = scaling(0.0, ScalingFailure::Warn); let result = solve_batch(warning_effects(), &ys, None, None, &cfg).expect("solve_batch"); assert!( @@ -76,3 +63,29 @@ fn free_solve_batch_surfaces_build_warnings() { let solver = Solver::new(warning_effects(), None, &cfg).expect("solver"); assert_eq!(result.warnings.as_slice(), solver.warnings()); } + +fn scaling(tolerance: f64, on_failure: ScalingFailure) -> PreconditionerConfig { + PreconditionerConfig::Additive { + local_solver: LocalSolverConfig { + scaling: ScalingConfig { + tolerance, + max_iterations: 0, + on_failure, + }, + ..LocalSolverConfig::default() + }, + } +} + +// The fixture's unscaled violation is about 0.078. +#[rstest::rstest] +#[case::within_tolerance(0.1, true)] +#[case::past_tolerance(0.05, false)] +fn only_a_violation_past_tolerance_fails_the_build(#[case] tolerance: f64, #[case] builds: bool) { + let built = Solver::new( + warning_effects(), + None, + scaling(tolerance, ScalingFailure::Error), + ); + assert_eq!(built.is_ok(), builds); +} diff --git a/python/within/README.md b/python/within/README.md index a323dc90..f004e254 100644 --- a/python/within/README.md +++ b/python/within/README.md @@ -112,7 +112,7 @@ seconds as `.build_duration_seconds`; both are preserved across pickling and reu | `LocalSolverConfig(approx_chol?, schur?, dense_threshold?, scaling?)` | Schur complement reduction with approximate Cholesky on the reduced system. Default local solver. Import from `within.config`. Omit `schur` for the library default (approximate); pass `Schur.exact()` for the exact complement. | | `Schur.approximate(config?)` / `Schur.exact()` | Schur-reduction mode for `LocalSolverConfig`. | | `ApproxCholConfig(seed=0, split_merge=None)` | Approximate Cholesky parameters. | -| `ApproxSchurConfig(seed=0, split=1)` | Approximate Schur complement sampling parameters. | +| `ApproxSchurConfig(split=1)` | Approximate Schur complement sampling parameters; draws from `ApproxCholConfig.seed`. | ### Result types diff --git a/python/within/_within.pyi b/python/within/_within.pyi index 722db540..482c5b3d 100644 --- a/python/within/_within.pyi +++ b/python/within/_within.pyi @@ -428,13 +428,14 @@ class ApproxCholConfig: def __init__(self, seed: int = 0, split_merge: int | None = None) -> None: ... class ApproxSchurConfig: - """Configuration for approximate Schur complement via clique-tree sampling.""" + """Configuration for approximate Schur complement via clique-tree sampling. + + The sampler draws from ``ApproxCholConfig.seed``. + """ - @property - def seed(self) -> int: ... @property def split(self) -> int: ... - def __init__(self, seed: int = 0, split: int = 1) -> None: ... + def __init__(self, split: int = 1) -> None: ... class Schur: """Schur-complement reduction mode for :class:`LocalSolverConfig`. diff --git a/tests/test_preconditioner.py b/tests/test_preconditioner.py index 2bf42f01..f563b316 100644 --- a/tests/test_preconditioner.py +++ b/tests/test_preconditioner.py @@ -84,7 +84,6 @@ def test_approx_chol_config_custom(self): def test_approx_schur_config_defaults(self): cfg = ApproxSchurConfig() - assert cfg.seed == 0 assert cfg.split == 1 def test_schur_complement_defaults(self): @@ -92,7 +91,7 @@ def test_schur_complement_defaults(self): assert sc.dense_threshold == 24 assert sc.approx_chol.seed == 0 assert sc.approx_chol.split_merge == 2 - assert repr(sc.schur) == "Schur.approximate(seed=0, split=1)" + assert repr(sc.schur) == "Schur.approximate(split=1)" assert sc.scaling.on_failure == "warn" def test_schur_mode_spellings(self): @@ -183,7 +182,7 @@ def test_preconditioner_exposes_tuned_config(self, problem): requested = PreconditionerConfig.Additive( local_solver=LocalSolverConfig( approx_chol=ApproxCholConfig(seed=41, split_merge=3), - schur=Schur.approximate(ApproxSchurConfig(seed=17, split=2)), + schur=Schur.approximate(ApproxSchurConfig(split=2)), dense_threshold=0, ), ) @@ -197,7 +196,7 @@ def test_preconditioner_exposes_tuned_config(self, problem): assert isinstance(actual, PreconditionerConfig.Additive) assert actual.local_solver.approx_chol.seed == 41 assert actual.local_solver.approx_chol.split_merge == 3 - assert repr(actual.local_solver.schur) == "Schur.approximate(seed=17, split=2)" + assert repr(actual.local_solver.schur) == "Schur.approximate(split=2)" assert actual.local_solver.dense_threshold == 0 assert actual.local_solver.scaling.on_failure == "warn"