From 7068444e4ba1c32781c60fe8d1598eca6d4109c8 Mon Sep 17 00:00:00 2001 From: Kristof Schroeder Date: Fri, 30 Oct 2026 17:15:49 +0100 Subject: [PATCH 1/3] feat!: hand each channel pair's factor-width-2 V to approx-chol 0.7 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Each channel pair's Gram matrix is V Vᵀ, one column per observation with its two levels' loadings scaled by √w. approx-chol 0.7 takes that bipartite factor-width-2 input and does the SGDD → SDDM reduction itself, so each of its irreducible blocks becomes a Schwarz subdomain behind BlockSolver. ApproxSchurConfig::seed is removed: the sampler draws from ApproxCholConfig::seed. The v21 wire fixture is regenerated. Refs #481 --- CHANGELOG.md | 6 +- Cargo.lock | 5 +- README.md | 4 +- benchmarks/_framework.py | 4 +- benchmarks/suites/ac_comparison.py | 2 +- benchmarks/suites/fixest_comparison.py | 4 +- crates/within-py/src/config.rs | 15 +- crates/within/Cargo.toml | 2 +- crates/within/src/block_solver.rs | 60 +++++ crates/within/src/config.rs | 62 ++--- crates/within/src/domain.rs | 7 +- crates/within/src/domain/factor_pairs.rs | 249 ++++++------------ crates/within/src/domain/prepared.rs | 9 +- crates/within/src/error.rs | 4 +- crates/within/src/lib.rs | 3 +- crates/within/src/linalg.rs | 8 - crates/within/src/operator/schwarz.rs | 15 +- crates/within/src/operator/schwarz/tests.rs | 129 +++------ crates/within/src/solver/tests.rs | 25 +- .../fixtures/preconditioner_v21.postcard | Bin 498 -> 386 bytes crates/within/tests/solver.rs | 2 +- crates/within/tests/warnings.rs | 6 +- python/within/README.md | 2 +- python/within/_within.pyi | 9 +- tests/test_preconditioner.py | 7 +- 25 files changed, 257 insertions(+), 382 deletions(-) create mode 100644 crates/within/src/block_solver.rs 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_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..fc8c1010 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 block dimension factored by exact dense Cholesky; `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/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/factor_pairs.rs b/crates/within/src/domain/factor_pairs.rs index b96e0f80..343c79af 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,90 @@ 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) + { + 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 +291,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/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/fixtures/preconditioner_v21.postcard b/crates/within/tests/fixtures/preconditioner_v21.postcard index ee2f4bb5f1bbef5d2965e19a143f625bb3b425f5..be7e672fd1f47f029edb574ace9e11bfa700856a 100644 GIT binary patch literal 386 zcmZQzVr5}vV_;-p>}=@iVqsuqVrF4rWMpP!VPIlm14}S4vN3?cgZ;$POw3S~_Lwvi z(5_H0cxul8vK*)w#d2oW#6>yRrdB+$pJVg>%b5qa?O}X~VIQD|ey~S%E z?1#`-_Ol!*GT8X{u0527iQ7Zf*)uVl3U?YN_I$SAm@VdWP=wV1#%E+?VPW9<_a6+P zwlOfk_~8)8?}vEe=zhVcrc<38Q%i literal 498 zcmZQzVr5}vV_;-pVPIlm12KS*jR6ck*b_@LGchwTgLE?i5y`rln3*9KJ%Cz>M7K2b zbg=*pWn^U3Q{)usoADW>2IOX-l_+lR^Z3~OdfS8j5c79m&J)i9xRyq0{6k&CM@fo3RW@TjLI%{2@HgCZ* z`?O#G?$&DG+yB_cMt#CvnEx0U_&hz5gl&)PSJ`gm`iU2&Zr}d&#Dp7?e|GJ^zQ^LY zw!?${N2b*V@?N@Y&%ngU>J^)FZED36`_d(szMOe*+aAV;M%I29%?^qV21bdg*VJBc P>j~O52rw`_oGb_c{fe7U 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..ff54eb08 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![ 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" From e6bd37e1be8899547fa59a2f7743cfe67c35acf3 Mon Sep 17 00:00:00 2001 From: Kristof Schroeder Date: Fri, 30 Oct 2026 17:17:11 +0100 Subject: [PATCH 2/3] refactor: delete within's own SGDD reduction Block elimination, the cross tables, the SDDM conversion with its scaling and grounding, the Schur assembly and the CSR block are unreferenced since the reduction moved into approx-chol. Refs #481 --- crates/within/src/block_elim.rs | 17 - crates/within/src/block_elim/csr_matrix.rs | 82 --- crates/within/src/block_elim/factor.rs | 248 --------- crates/within/src/block_elim/schur.rs | 473 ------------------ crates/within/src/block_elim/schur/tests.rs | 205 -------- crates/within/src/block_elim/solver.rs | 468 ----------------- crates/within/src/block_elim/solver/tests.rs | 442 ---------------- crates/within/src/block_elim/tests.rs | 34 -- crates/within/src/csr_block.rs | 368 -------------- crates/within/src/domain/cross_tab.rs | 237 --------- .../within/src/domain/cross_tab/accumulate.rs | 270 ---------- crates/within/src/domain/cross_tab/tests.rs | 277 ---------- .../src/domain/factor_pairs/grounding.rs | 53 -- crates/within/src/domain/factor_pairs/sddm.rs | 383 -------------- .../src/domain/factor_pairs/sddm/scaling.rs | 401 --------------- .../src/domain/factor_pairs/sddm/tests.rs | 340 ------------- 16 files changed, 4298 deletions(-) delete mode 100644 crates/within/src/block_elim.rs delete mode 100644 crates/within/src/block_elim/csr_matrix.rs delete mode 100644 crates/within/src/block_elim/factor.rs delete mode 100644 crates/within/src/block_elim/schur.rs delete mode 100644 crates/within/src/block_elim/schur/tests.rs delete mode 100644 crates/within/src/block_elim/solver.rs delete mode 100644 crates/within/src/block_elim/solver/tests.rs delete mode 100644 crates/within/src/block_elim/tests.rs delete mode 100644 crates/within/src/csr_block.rs delete mode 100644 crates/within/src/domain/cross_tab.rs delete mode 100644 crates/within/src/domain/cross_tab/accumulate.rs delete mode 100644 crates/within/src/domain/cross_tab/tests.rs delete mode 100644 crates/within/src/domain/factor_pairs/grounding.rs delete mode 100644 crates/within/src/domain/factor_pairs/sddm.rs delete mode 100644 crates/within/src/domain/factor_pairs/sddm/scaling.rs delete mode 100644 crates/within/src/domain/factor_pairs/sddm/tests.rs 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/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/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/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); - } - } - } -} From 575917903ba37af1bbaca9cd5863cd6b80e40e23 Mon Sep 17 00:00:00 2001 From: Kristof Schroeder Date: Sat, 31 Oct 2026 17:05:07 +0100 Subject: [PATCH 3/3] fix: fail a scaling only past its tolerance; test the tolerance < 1 boundary approx-chol 0.7 reports any deficit past rounding, so filter on violation > tolerance as before. Document dense_threshold as the Cholesky dimension with a signed block's cover sheets counted. --- crates/within/src/config.rs | 2 +- crates/within/src/domain/factor_pairs.rs | 1 + crates/within/tests/error_paths.rs | 2 +- crates/within/tests/warnings.rs | 43 +++++++++++++++--------- 4 files changed, 31 insertions(+), 17 deletions(-) diff --git a/crates/within/src/config.rs b/crates/within/src/config.rs index fc8c1010..a2684b79 100644 --- a/crates/within/src/config.rs +++ b/crates/within/src/config.rs @@ -47,7 +47,7 @@ pub struct LocalSolverConfig { pub approx_chol: ApproxCholConfig, /// Schur-complement reduction mode (default: approximate). pub schur: SchurMode, - /// Largest block dimension factored by exact dense Cholesky; `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, diff --git a/crates/within/src/domain/factor_pairs.rs b/crates/within/src/domain/factor_pairs.rs index 343c79af..ec55e8cd 100644 --- a/crates/within/src/domain/factor_pairs.rs +++ b/crates/within/src/domain/factor_pairs.rs @@ -109,6 +109,7 @@ fn factor_pair( .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 { 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/warnings.rs b/crates/within/tests/warnings.rs index ff54eb08..aadcc287 100644 --- a/crates/within/tests/warnings.rs +++ b/crates/within/tests/warnings.rs @@ -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); +}