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