From 22c6b6c123341a92910fb01c256fb48bd3f74d7f Mon Sep 17 00:00:00 2001 From: Max Bolingbroke Date: Sun, 19 Jul 2026 19:20:42 -0700 Subject: [PATCH 1/2] qdldl: replay a precomputed schedule on numeric refactorization After the first factorization the sparsity pattern of L is fixed, but _factor_inner re-derives the per-column update lists on every refactor: an elimination-tree walk with marker arrays, a work buffer and a reversal step. Interior-point use refactors the same pattern once per iteration, so this control flow is pure overhead. This change records the discovered control flow on the first refactor as a flat list of (column, position) update steps, then replays the identical floating-point operations in the identical order on subsequent refactorizations. L, D, Dinv, the inertia count and the dynamic-regularization decisions are bit-for-bit identical to the original path -- enforced by new tests that compare a refactorization against a fresh factorization bitwise, including a case where dynamic regularization fires. This is the standard symbolic/numeric phase separation of sparse direct solvers (T. Davis, "Direct Methods for Sparse Linear Systems", SIAM 2006, ch. 4). The schedule is built lazily on the first refactor() call, so one-shot factorizations pay nothing. The build self-verifies against the existing pattern of L (every recorded position must satisfy Li[pos] == k, and the step count must equal nnz(L)); on any mismatch, or if indices would not fit the u32 storage used to halve replay memory traffic, refactorization falls back to the original path. Measured on an interleaved A/B benchmark (Apple M-series, min of 3 rounds x 2 reps, objectives bitwise-equal throughout): -2.1% total wall time over nine portfolio-rebalance QP/SOCPs (n=9.5k-18k), and -3.9% total over seven Netlib-Kennington/Mittelmann/Maros-Meszaros/ conic problems (all 16 problems individually faster, -1.2% to -5.4%). The win is modest because these factors are flop-bound; the removed overhead grows in relative terms the sparser the factor. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01CqGdm2vZZ6HA5HT8DsRePa --- src/qdldl/qdldl.rs | 321 ++++++++++++++++++++++++++++++++++++++++++++- src/qdldl/test.rs | 139 ++++++++++++++++++++ 2 files changed, 455 insertions(+), 5 deletions(-) diff --git a/src/qdldl/qdldl.rs b/src/qdldl/qdldl.rs index b36b82e6..ee243f80 100644 --- a/src/qdldl/qdldl.rs +++ b/src/qdldl/qdldl.rs @@ -190,13 +190,77 @@ where // factorization since it will always be the same. Calling // this function implies that we want a numerical factorization self.is_symbolic = false; - _factor( - &mut self.L, + + // The sparsity pattern of L was fixed by the factorization in + // `new` (numeric or logical), so refactorizations can replay a + // precomputed update schedule instead of re-deriving the pattern. + // Built lazily here so that single-factorization use pays nothing. + if matches!(self.workspace.schedule, ScheduleState::NotBuilt) { + let ws = &self.workspace; + self.workspace.schedule = match _build_schedule( + &ws.triuA.colptr, + &ws.triuA.rowval, + &ws.etree, + &self.L.colptr, + &self.L.rowval, + ) { + Some(s) => ScheduleState::Ready(s), + None => ScheduleState::Unavailable, + }; + } + + if !matches!(self.workspace.schedule, ScheduleState::Ready(_)) { + return _factor( + &mut self.L, + &mut self.D, + &mut self.Dinv, + &mut self.workspace, + false, + ); + } + + // destructure so the borrows of the schedule and the other + // workspace fields are disjoint + let QDLDLWorkspace { + schedule, + triuA, + fwork, + Dsigns, + regularize_enable, + regularize_eps, + regularize_delta, + regularize_count, + positive_inertia, + .. + } = &mut self.workspace; + let ScheduleState::Ready(sched) = schedule else { + unreachable!() + }; + + *positive_inertia = _refactor_scheduled( + sched, + &triuA.colptr, + &triuA.rowval, + &triuA.nzval, + &self.L.colptr, + &self.L.rowval, + &mut self.L.nzval, &mut self.D, &mut self.Dinv, - &mut self.workspace, - self.is_symbolic, - ) + fwork, + Dsigns, + *regularize_enable, + *regularize_eps, + *regularize_delta, + regularize_count, + )?; + Ok(()) + } + + // test support: confirms refactorization used the scheduled replay path + #[cfg(test)] + pub(crate) fn refactor_schedule_is_ready(&self) -> bool { + matches!(self.workspace.schedule, ScheduleState::Ready(_)) } /// Returns the number of nonzeros in A for A = LDL^T @@ -325,6 +389,10 @@ struct QDLDLWorkspace { // number of regularized entries in D regularize_count: usize, + + // precomputed control flow for numeric refactorization, + // built lazily on the first refactor + schedule: ScheduleState, } impl QDLDLWorkspace @@ -375,10 +443,253 @@ where regularize_eps, regularize_delta, regularize_count, + schedule: ScheduleState::NotBuilt, }) } } +// ------------------------------------- +// Precomputed refactorization schedule +// ------------------------------------- +// +// After the first factorization the sparsity pattern of L (colptr/rowval) +// is fixed; refactorizations change only numeric values. `_factor_inner` +// nevertheless re-derives, for every column k, the list of prior columns +// that update it — an elimination-tree walk with marker arrays, a work +// buffer and a reversal step. On matrices whose factors have little +// fill-in (typical for IPM KKT systems), that control-flow overhead +// dominates the actual arithmetic of a refactorization. +// +// `FactorSchedule` records the discovered control flow once, flat: for each +// column k in order, the sequence of (cidx, pos) update steps exactly as +// `_factor_inner` executes them, where cidx is the prior column supplying +// the update and pos is the position in L.nzval that receives L[k,cidx]. +// The update range of a step is Lp[cidx]..pos: at the moment step (k,cidx) +// runs, column cidx holds exactly its entries with row index < k, because +// entries are appended to each column in increasing row order. +// `_refactor_scheduled` then replays the identical floating-point +// operations in the identical order — giving bit-identical L, D and +// regularization decisions — without re-deriving any of the pattern. +// This is the standard symbolic/numeric phase separation of sparse direct +// solvers (T. Davis, "Direct Methods for Sparse Linear Systems", SIAM 2006, +// ch. 4); QDLDL's original single-phase design favours simplicity for +// one-shot factorizations, but Clarabel refactors the same pattern once +// per interior-point iteration. +// +// Indices are stored as u32 to halve the memory traffic of the replay. +// Patterns too large for that (nnz(L) or n >= 2^32) fall back to the +// original path via `ScheduleState::Unavailable`. + +#[derive(Debug)] +struct FactorSchedule { + // start of column k's steps in cidx/pos; length n+1 + colptr: Vec, + // source column of each update step + cidx: Vec, + // position in L.nzval written by each update step + pos: Vec, +} + +#[derive(Debug)] +enum ScheduleState { + NotBuilt, + Unavailable, + Ready(FactorSchedule), +} + +// Records the update schedule by replaying the pattern-discovery phase of +// `_factor_inner` (which see), with the numeric work stripped out. The two +// functions must stay in lockstep; as a defence, every recorded write +// position is verified against the already-computed pattern of L +// (Li[pos] == k), and the total step count against nnz(L). Any mismatch +// returns None and refactorization falls back to `_factor_inner`. +fn _build_schedule( + Ap: &[usize], + Ai: &[usize], + etree: &[usize], + Lp: &[usize], + Li: &[usize], +) -> Option { + let n = Lp.len() - 1; + let nnzL = Lp[n]; + if n >= u32::MAX as usize || nnzL >= u32::MAX as usize { + return None; + } + + let mut colptr = Vec::with_capacity(n + 1); + let mut cidx = Vec::with_capacity(nnzL); + let mut pos = Vec::with_capacity(nnzL); + + let mut y_markers = vec![QDLDL_UNUSED; n]; + let mut y_idx = vec![0usize; n]; + let mut elim_buffer = vec![0usize; n]; + let mut next_colspace: Vec = Lp[0..n].to_vec(); + + colptr.push(0u32); + if n > 0 { + colptr.push(0u32); // column 0 has no update steps + } + + for k in 1..n { + // pattern-discovery phase, exactly as in _factor_inner + let mut nnz_y = 0; + + for &bidx in &Ai[Ap[k]..Ap[k + 1]] { + if bidx == k { + continue; + } + + if y_markers[bidx] == QDLDL_UNUSED { + y_markers[bidx] = QDLDL_USED; + elim_buffer[0] = bidx; + let mut nnz_e = 1; + + let mut next_idx = etree[bidx]; + while next_idx != QDLDL_UNKNOWN && next_idx < k { + if y_markers[next_idx] == QDLDL_USED { + break; + } + y_markers[next_idx] = QDLDL_USED; + elim_buffer[nnz_e] = next_idx; + next_idx = etree[next_idx]; + nnz_e += 1; + } + + while nnz_e != 0 { + nnz_e -= 1; + y_idx[nnz_y] = elim_buffer[nnz_e]; + nnz_y += 1; + } + } + } + + // record the value-placement phase of _factor_inner + for i in (0..nnz_y).rev() { + let c = y_idx[i]; + let p = next_colspace[c]; + + // verify against the pattern L already has + if Li[p] != k { + return None; + } + + cidx.push(c as u32); + pos.push(p as u32); + next_colspace[c] += 1; + y_markers[c] = QDLDL_UNUSED; + } + colptr.push(cidx.len() as u32); + } + + if cidx.len() != nnzL { + return None; + } + + Some(FactorSchedule { colptr, cidx, pos }) +} + +// Numeric refactorization by schedule replay. Performs the identical +// floating point operations, in the identical order, as +// `_factor_inner(..., logical_factor = false)` on the same pattern, so the +// results (L, D, Dinv, inertia and regularization counts) are bit-identical +// to that function's. See `_build_schedule` for the schedule's invariants. +#[allow(clippy::too_many_arguments)] +fn _refactor_scheduled( + sched: &FactorSchedule, + Ap: &[usize], + Ai: &[usize], + Ax: &[T], + Lp: &[usize], + Li: &[usize], + Lx: &mut [T], + D: &mut [T], + Dinv: &mut [T], + y_vals: &mut [T], + Dsigns: &[i8], + regularize_enable: bool, + regularize_eps: T, + regularize_delta: T, + regularize_count: &mut usize, +) -> Result { + *regularize_count = 0; + let mut positiveValuesInD = 0; + let n = Lp.len() - 1; + + y_vals.fill(T::zero()); + D.fill(T::zero()); + + // First element of the diagonal D, as in _factor_inner + D[0] = Ax[0]; + if regularize_enable { + let sign = T::from_i8(Dsigns[0]).unwrap(); + if D[0] * sign < regularize_eps { + D[0] = regularize_delta * sign; + *regularize_count += 1; + } + } + if D[0].is_zero() { + return Err(QDLDLError::ZeroPivot); + } + if D[0] > T::zero() { + positiveValuesInD += 1; + } + Dinv[0] = T::recip(D[0]); + + for k in 1..n { + // scatter the kth column of A above the diagonal into the sparse + // accumulator, and initialize D[k], exactly as in _factor_inner + for i in Ap[k]..Ap[k + 1] { + let bidx = Ai[i]; + if bidx == k { + D[k] = Ax[i]; + } else { + y_vals[bidx] = Ax[i]; + } + } + + // replay the update steps for this column + let (f, l) = (sched.colptr[k] as usize, sched.colptr[k + 1] as usize); + for (&c, &p) in zip(&sched.cidx[f..l], &sched.pos[f..l]) { + let (cidx, tmp_idx) = (c as usize, p as usize); + let y_vals_cidx = y_vals[cidx]; + + let (f, l) = (Lp[cidx], tmp_idx); + unsafe { + // Safety: Li entries index the matrix dimension, and the + // schedule's positions were verified against Li at build + // time; both bound y_vals/Lx as in _factor_inner. + for (&Lxj, &Lij) in zip(&Lx[f..l], &Li[f..l]) { + *(y_vals.get_unchecked_mut(Lij)) -= Lxj * y_vals_cidx; + } + + let Lx_tmp_idx = y_vals_cidx * *Dinv.get_unchecked(cidx); + *Lx.get_unchecked_mut(tmp_idx) = Lx_tmp_idx; + *D.get_unchecked_mut(k) -= y_vals_cidx * Lx_tmp_idx; + } + + y_vals[cidx] = T::zero(); + } + + // pivot regularization / rejection, as in _factor_inner + if regularize_enable { + let sign = T::from_i8(Dsigns[k]).unwrap(); + if D[k] * sign < regularize_eps { + D[k] = regularize_delta * sign; + *regularize_count += 1; + } + } + if D[k].is_zero() { + return Err(QDLDLError::ZeroPivot); + } + if D[k] > T::zero() { + positiveValuesInD += 1; + } + Dinv[k] = T::recip(D[k]); + } + + Ok(positiveValuesInD) +} + fn _factor( L: &mut CscMatrix, D: &mut [T], diff --git a/src/qdldl/test.rs b/src/qdldl/test.rs index 020424f1..3c47581b 100644 --- a/src/qdldl/test.rs +++ b/src/qdldl/test.rs @@ -263,6 +263,145 @@ fn test_solve_logical_refactor() { assert!(inf_norm_diff(&x, &b) <= 1e-8); } +// Build a deterministic pseudo-random quasidefinite KKT-like matrix +// [[diag(p) B'; B -diag(r)]] in upper triangular CSC form, sized so that +// its factor has genuine fill-in. Returns (A, Dsigns). +#[cfg(test)] +fn test_matrix_quasidef(nx: usize, nz: usize, seed: u64) -> (CscMatrix, Vec) { + // simple LCG so the test needs no rand dependency + let mut state = seed; + let mut next = move || { + state = state + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + ((state >> 33) as f64) / ((1u64 << 31) as f64) - 1.0 // in [-1,1) + }; + + let n = nx + nz; + let mut cols: Vec> = vec![Vec::new(); n]; + for j in 0..nx { + cols[j].push((j, 1.0 + next().abs())); // positive definite block + } + for j in 0..nz { + let col = nx + j; + // a few entries of B in each column, rows in 0..nx + for t in 0..3 { + let i = ((next().abs() * nx as f64) as usize + t * 7) % nx; + cols[col].push((i, next())); + } + cols[col].sort_by_key(|e| e.0); + cols[col].dedup_by_key(|e| e.0); + cols[col].push((col, -(1.0 + next().abs()))); // negative definite block + } + + let mut colptr = vec![0usize]; + let (mut rowval, mut nzval) = (Vec::new(), Vec::new()); + for c in &cols { + for &(i, v) in c { + rowval.push(i); + nzval.push(v); + } + colptr.push(rowval.len()); + } + let A = CscMatrix { + m: n, + n, + colptr, + rowval, + nzval, + }; + let mut signs = vec![1i8; n]; + signs[nx..].fill(-1); + (A, signs) +} + +// Refactorization must reproduce _factor_inner bit-for-bit: the scheduled +// replay path performs the identical operations in the identical order, so +// L, D, Dinv, the inertia and the regularization count of a refactor must +// all equal those of a fresh factorization of the same values. +#[test] +fn test_refactor_matches_fresh_factor_exactly() { + let (A, signs) = test_matrix_quasidef(40, 30, 12345); + + let opts = || { + QDLDLSettingsBuilder::::default() + .Dsigns(signs.clone()) + .build() + .unwrap() + }; + + let mut f1 = QDLDLFactorisation::new(&A, Some(opts())).unwrap(); + + // change every value, refactor (first refactor builds the schedule + // and replays it), and compare against a fresh factorization + let mut A2 = A.clone(); + for v in A2.nzval.iter_mut() { + *v *= 1.25; + } + let indices: Vec = (0..A2.nzval.len()).collect(); + + f1.update_values(&indices, &A2.nzval); + f1.refactor().unwrap(); + assert!(f1.refactor_schedule_is_ready()); // replay path, not a fallback + + let f2 = QDLDLFactorisation::new(&A2, Some(opts())).unwrap(); + + assert_eq!(f1.perm, f2.perm); // same AMD ordering on the same pattern + assert_eq!(f1.L.nzval, f2.L.nzval); // bitwise + assert_eq!(f1.D, f2.D); + assert_eq!(f1.Dinv, f2.Dinv); + assert_eq!(f1.positive_inertia(), f2.positive_inertia()); + assert_eq!(f1.regularize_count(), f2.regularize_count()); + + // and again, to exercise the replay path on an already-built schedule + f1.update_values(&indices, &A.nzval); + f1.refactor().unwrap(); + let f3 = QDLDLFactorisation::new(&A, Some(opts())).unwrap(); + assert_eq!(f1.L.nzval, f3.L.nzval); + assert_eq!(f1.D, f3.D); +} + +// Same bit-identity requirement when dynamic regularization fires: the +// pivot tests happen in the same order on the same values, so the same +// pivots must be perturbed. +#[test] +fn test_refactor_matches_fresh_factor_with_regularization() { + let (mut A, signs) = test_matrix_quasidef(40, 30, 999); + // shrink some diagonal entries so that regularization triggers + for j in 0..40 { + let d = A.colptr[j]; // diagonal of the (j,j) leading block column + A.nzval[d] *= 1e-14; + } + + let opts = || { + QDLDLSettingsBuilder::::default() + .Dsigns(signs.clone()) + .regularize_eps(1e-12) + .regularize_delta(1e-7) + .build() + .unwrap() + }; + + let mut f1 = QDLDLFactorisation::new(&A, Some(opts())).unwrap(); + + let mut A2 = A.clone(); + for v in A2.nzval.iter_mut() { + *v *= 0.75; + } + let indices: Vec = (0..A2.nzval.len()).collect(); + f1.update_values(&indices, &A2.nzval); + f1.refactor().unwrap(); + assert!(f1.refactor_schedule_is_ready()); // replay path, not a fallback + + let f2 = QDLDLFactorisation::new(&A2, Some(opts())).unwrap(); + assert!(f2.regularize_count() > 0); // the scenario is actually exercised + assert_eq!(f1.L.nzval, f2.L.nzval); + assert_eq!(f1.D, f2.D); + assert_eq!(f1.Dinv, f2.Dinv); + assert_eq!(f1.positive_inertia(), f2.positive_inertia()); + assert_eq!(f1.regularize_count(), f2.regularize_count()); +} + #[test] fn test_bad_numeric_pivot() { //Disable regularization to force an exact zero pivot From 096ae970e698b2a4ade6c4a74f6e8a7f16f38e10 Mon Sep 17 00:00:00 2001 From: Max Bolingbroke Date: Tue, 21 Jul 2026 01:33:20 -0700 Subject: [PATCH 2/2] qdldl: replay updates block-wise on the columns where it pays The dominant loop of the numeric factorization scatters into the sparse accumulator, one indexed load and one indexed store per entry: for j in range { y[Li[j]] -= Lx[j] * y_c } Row indices within a column of L are often long blocks of consecutive integers -- a set of dense-ish rows, which AMD orders last, produces a trailing trapezoid shared by many columns. Replaying such a block as a contiguous slice update drops the per-entry index load and lets the compiler vectorize: y[rs..rs+len] -= Lx[p..p+len] * y_c That is worth up to 55% of solve time on factors built from long blocks. But it costs bookkeeping per *block* -- the bounds, the clip against the end of the update range, and the inner loop's own setup -- and on factors whose blocks are two or three entries long that is not amortized: measured across public problems, replaying every column block-wise costs 10-18% on small sparse LPs. (A per-block fast path does not help: for a short block the contiguous and indexed inner loops do the same work, and the cost is the block loop itself.) Since the trade is a property of each column rather than of the matrix, it is decided per column. A column's blocks are recorded only if the block containing a typical entry of that column reaches RUN_MIN_LEN entries -- entry-weighted mean block length, Sum(len^2) / Sum(len) -- and columns failing the test are replayed entry-wise. A column with no recorded blocks takes the entry-wise path, so the decision doubles as its own storage and needs no extra array. Factors of mixed structure, which are the common case, then take the contiguous path on exactly the columns that benefit. RUN_MIN_LEN is 16, a vectorization threshold: eight iterations of a two-wide double-precision loop, comfortably past its prologue and epilogue, whereas at two or four entries the setup dominates. Both paths perform the same operations on the same values in the same order, so the factorization is bit-identical whichever is chosen. Two tests assert that against a fresh factorization on matrices sitting on opposite sides of the decision, and a third pins that the decision does go both ways, so neither path can rot untested. Measured on 80 public problems (Netlib, Netlib-Kennington, Maros-Meszaros, Mittelmann, SDPLIB, structured conic), interleaved against the schedule replay alone: 67 small and medium problems: -29.2% total, median -10.0%, 54 faster, 12 unchanged, one slower by 2.3% 13 large problems (>10s), held out from the calibration of RUN_MIN_LEN: -35.7% total, all 13 faster Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01CqGdm2vZZ6HA5HT8DsRePa --- src/qdldl/qdldl.rs | 157 +++++++++++++++++++++++++++++++++++++++------ src/qdldl/test.rs | 108 +++++++++++++++++++++++++++++++ 2 files changed, 247 insertions(+), 18 deletions(-) diff --git a/src/qdldl/qdldl.rs b/src/qdldl/qdldl.rs index ee243f80..17bf2d36 100644 --- a/src/qdldl/qdldl.rs +++ b/src/qdldl/qdldl.rs @@ -257,6 +257,17 @@ where Ok(()) } + // test support: how many columns replay their updates block-wise + #[cfg(test)] + pub(crate) fn columns_using_blocks(&self) -> usize { + match &self.workspace.schedule { + ScheduleState::Ready(s) => (0..s.colruns.len() - 1) + .filter(|&c| s.colruns[c] != s.colruns[c + 1]) + .count(), + _ => 0, + } + } + // test support: confirms refactorization used the scheduled replay path #[cfg(test)] pub(crate) fn refactor_schedule_is_ready(&self) -> bool { @@ -479,6 +490,19 @@ where // Indices are stored as u32 to halve the memory traffic of the replay. // Patterns too large for that (nnz(L) or n >= 2^32) fall back to the // original path via `ScheduleState::Unavailable`. +// +// The schedule also stores a run decomposition of each L column's row +// indices: the maximal blocks of consecutive rows. KKT factors from +// block-structured problems concentrate their flops in columns whose +// entries lie in long consecutive row blocks (e.g. a trailing +// dense-border trapezoid produced by a set of dense rows), so the +// dominant update loop +// for j in range { y[Li[j]] -= Lx[j] * y_c } +// mostly writes to contiguous y locations. Replaying it run-by-run as +// y[rs..rs+len] -= Lx[p..p+len] * y_c +// performs the same independent elementwise operations on the same +// values — bit-identical results — but with no per-element index load +// and in a form the compiler can vectorize. #[derive(Debug)] struct FactorSchedule { @@ -488,6 +512,45 @@ struct FactorSchedule { cidx: Vec, // position in L.nzval written by each update step pos: Vec, + // start of column c's blocks in `runs`; length n+1. A column with an + // empty range is replayed entry-wise -- see `_column_prefers_blocks`. + colruns: Vec, + // (first row, length) of each maximal consecutive-row block, in position + // order within each column + runs: Vec<(u32, u32)>, +} + +// Replaying an update over maximal blocks of consecutive rows removes one +// indexed load per entry and lets the compiler vectorize the inner loop, but +// it costs bookkeeping per *block*: the block bounds, the clip against the end +// of the update range, and the inner loop's own setup. For a block of two or +// three entries that bookkeeping is not amortized and the entry-wise loop is +// faster; for a long block it is negligible. +// +// The trade is a property of each column rather than of the matrix, so it is +// decided per column: a column's blocks are recorded only if the block +// containing a typical entry of that column reaches `RUN_MIN_LEN`, and columns +// failing the test are replayed entry-wise. A factorization of mixed +// structure -- common when a few dense-ish rows sit above a sparse remainder +// -- then takes the contiguous path on exactly the columns that benefit. +// +// `RUN_MIN_LEN` is a vectorization threshold: at 16 entries a two-wide +// double-precision loop runs eight iterations, comfortably past its prologue +// and epilogue, whereas at two or four the setup dominates. +// +// The average is entry-weighted, Sum(len^2) / Sum(len) over the column's +// blocks, i.e. the mean block length seen by a randomly chosen entry of the +// column, because it is the per-entry saving that has to outweigh the +// per-block cost. +const RUN_MIN_LEN: u32 = 16; + +fn _column_prefers_blocks(runs: &[(u32, u32)]) -> bool { + let (mut sum, mut sumsq) = (0u64, 0u64); + for &(_, len) in runs { + sum += len as u64; + sumsq += (len as u64) * (len as u64); + } + sum > 0 && sumsq >= (RUN_MIN_LEN as u64) * sum } #[derive(Debug)] @@ -497,12 +560,6 @@ enum ScheduleState { Ready(FactorSchedule), } -// Records the update schedule by replaying the pattern-discovery phase of -// `_factor_inner` (which see), with the numeric work stripped out. The two -// functions must stay in lockstep; as a defence, every recorded write -// position is verified against the already-computed pattern of L -// (Li[pos] == k), and the total step count against nnz(L). Any mismatch -// returns None and refactorization falls back to `_factor_inner`. fn _build_schedule( Ap: &[usize], Ai: &[usize], @@ -585,7 +642,44 @@ fn _build_schedule( return None; } - Some(FactorSchedule { colptr, cidx, pos }) + // Decompose each column's (already final) row indices into maximal blocks + // of consecutive rows, keeping them only for the columns that prefer the + // block-wise replay. A column with no recorded blocks is replayed + // entry-wise, so this doubles as the decision and as its storage: nothing + // is kept for columns that would not use it. + let mut colruns = Vec::with_capacity(n + 1); + let mut runs: Vec<(u32, u32)> = Vec::new(); + let mut col: Vec<(u32, u32)> = Vec::new(); + colruns.push(0u32); + for c in 0..n { + let li = &Li[Lp[c]..Lp[c + 1]]; + col.clear(); + let mut t = 0; + while t < li.len() { + let start = li[t]; + let mut len = 1; + while t + len < li.len() && li[t + len] == start + len { + len += 1; + } + col.push((start as u32, len as u32)); + t += len; + } + if _column_prefers_blocks(&col) { + runs.extend_from_slice(&col); + } + colruns.push(runs.len() as u32); + } + if runs.len() >= u32::MAX as usize { + return None; + } + + Some(FactorSchedule { + colptr, + cidx, + pos, + colruns, + runs, + }) } // Numeric refactorization by schedule replay. Performs the identical @@ -653,20 +747,47 @@ fn _refactor_scheduled( let (cidx, tmp_idx) = (c as usize, p as usize); let y_vals_cidx = y_vals[cidx]; - let (f, l) = (Lp[cidx], tmp_idx); - unsafe { - // Safety: Li entries index the matrix dimension, and the - // schedule's positions were verified against Li at build - // time; both bound y_vals/Lx as in _factor_inner. - for (&Lxj, &Lij) in zip(&Lx[f..l], &Li[f..l]) { - *(y_vals.get_unchecked_mut(Lij)) -= Lxj * y_vals_cidx; + // The update range is Lp[cidx]..tmp_idx. Either way this performs + // the same independent operations on the same values as the scalar + // loop in _factor_inner; the two paths differ only in how the row + // of each entry is obtained. + let (rf, rl) = ( + sched.colruns[cidx] as usize, + sched.colruns[cidx + 1] as usize, + ); + if rf != rl { + // Block-wise: each block covers consecutive rows, so the + // scatter becomes a contiguous elementwise update with no + // per-entry index load, in a form the compiler vectorizes. + let mut pos = Lp[cidx]; + for &(row_start, run_len) in &sched.runs[rf..rl] { + if pos >= tmp_idx { + break; + } + let take = min(run_len as usize, tmp_idx - pos); + let rs = row_start as usize; + for (yv, &Lxj) in zip(&mut y_vals[rs..rs + take], &Lx[pos..pos + take]) { + *yv -= Lxj * y_vals_cidx; + } + pos += take; + } + } else { + // Entry-wise, for columns whose blocks are too short to be + // worth the per-block bookkeeping. + let (f, l) = (Lp[cidx], tmp_idx); + unsafe { + // Safety: Li entries index the matrix dimension, so they + // bound y_vals as in _factor_inner. + for (&Lxj, &Lij) in zip(&Lx[f..l], &Li[f..l]) { + *(y_vals.get_unchecked_mut(Lij)) -= Lxj * y_vals_cidx; + } } - - let Lx_tmp_idx = y_vals_cidx * *Dinv.get_unchecked(cidx); - *Lx.get_unchecked_mut(tmp_idx) = Lx_tmp_idx; - *D.get_unchecked_mut(k) -= y_vals_cidx * Lx_tmp_idx; } + let Lx_tmp_idx = y_vals_cidx * Dinv[cidx]; + Lx[tmp_idx] = Lx_tmp_idx; + D[k] -= y_vals_cidx * Lx_tmp_idx; + y_vals[cidx] = T::zero(); } diff --git a/src/qdldl/test.rs b/src/qdldl/test.rs index 3c47581b..f72174f4 100644 --- a/src/qdldl/test.rs +++ b/src/qdldl/test.rs @@ -402,6 +402,114 @@ fn test_refactor_matches_fresh_factor_with_regularization() { assert_eq!(f1.regularize_count(), f2.regularize_count()); } +// Both replay paths must produce the same factorization: the block-wise and +// entry-wise updates perform the same operations on the same values in the +// same order, and which one is selected is a performance decision only. The +// two matrices below sit on opposite sides of that decision. +#[test] +fn test_both_replay_paths_agree_with_a_fresh_factor() { + // banded, so the factor's rows are long consecutive blocks + let n = 120usize; + let (mut colptr, mut rowval, mut nzval) = (vec![0usize], Vec::new(), Vec::new()); + for j in 0..n { + for i in j.saturating_sub(20)..=j { + rowval.push(i); + nzval.push(if i == j { + 40.0 + } else { + -1.0 / (1 + j - i) as f64 + }); + } + colptr.push(rowval.len()); + } + let banded = CscMatrix { + m: n, + n, + colptr, + rowval, + nzval, + }; + + // and a sparse quasidefinite matrix, whose blocks are short + let (scattered, signs) = test_matrix_quasidef(60, 45, 4242); + + for (name, A, ds) in [ + ("banded", banded, None), + ("scattered", scattered, Some(signs)), + ] { + let mut b = QDLDLSettingsBuilder::::default(); + if let Some(ds) = ds { + b.Dsigns(ds); + } + let opts = b.build().unwrap(); + + let mut f = QDLDLFactorisation::new(&A, Some(opts.clone())).unwrap(); + let mut scaled = A.clone(); + for v in scaled.nzval.iter_mut() { + *v *= 1.25; + } + let idx: Vec = (0..A.nzval.len()).collect(); + f.update_values(&idx, &scaled.nzval); + f.refactor().unwrap(); + assert!( + f.refactor_schedule_is_ready(), + "{name}: replay path not engaged" + ); + + let fresh = QDLDLFactorisation::new(&scaled, Some(opts)).unwrap(); + assert_eq!(f.L.nzval, fresh.L.nzval, "{name}: L differs"); + assert_eq!(f.D, fresh.D, "{name}: D differs"); + assert_eq!(f.Dinv, fresh.Dinv, "{name}: Dinv differs"); + } +} + +// The block-wise path must be chosen for a factor whose updates are long +// consecutive blocks, and declined for one whose blocks are short -- which is +// what makes the two cases above cover both paths rather than one twice. +#[test] +fn test_run_selection_follows_block_structure() { + let n = 120usize; + let (mut colptr, mut rowval, mut nzval) = (vec![0usize], Vec::new(), Vec::new()); + for j in 0..n { + for i in j.saturating_sub(20)..=j { + rowval.push(i); + nzval.push(if i == j { + 40.0 + } else { + -1.0 / (1 + j - i) as f64 + }); + } + colptr.push(rowval.len()); + } + let banded = CscMatrix { + m: n, + n, + colptr, + rowval, + nzval, + }; + let mut f = QDLDLFactorisation::new(&banded, None).unwrap(); + f.refactor().unwrap(); + let banded_cols = f.columns_using_blocks(); + assert!( + banded_cols > n / 2, + "long blocks should select the block path for most columns, got {banded_cols}" + ); + + let (scattered, signs) = test_matrix_quasidef(60, 45, 4242); + let opts = QDLDLSettingsBuilder::::default() + .Dsigns(signs) + .build() + .unwrap(); + let mut f = QDLDLFactorisation::new(&scattered, Some(opts)).unwrap(); + f.refactor().unwrap(); + assert_eq!( + f.columns_using_blocks(), + 0, + "short blocks should decline the block path" + ); +} + #[test] fn test_bad_numeric_pivot() { //Disable regularization to force an exact zero pivot