svd: dqds takes its work arrays from the caller's MemStack (#156)

PATCHSET.md

@@ -106,6 +106,7 @@ gate has changed and two kinds of adoption are now in the ledger.
 | dqds for values-only singular values | LAPACK's `dgesdd(jobz='N')` never applies a rotation: `DBDSDC(compq='N')` -> `DLASDQ` -> `DBDSQR`'s `ROTATE`-false branch -> `DLASQ1`, the dqds iteration. The fork had no such path. `bidiag_svd::dqds` now applies the shifted qd transform to the squared factors (one division and one multiplication per entry per sweep) with scaling before squaring, a positivity-checked shift retry, splitting at zero off-diagonals and Li's deflation; a rank-deficient bidiagonal keeps the rotation path. Verified against `thin_svd` for n = 2..40 over four spectra at 1e-10. It is serial, and the divide-and-conquer it replaces threads its merge, so the route is by (budget, size): dqds at one thread or below n = 256, divide and conquer above. Campaign `svd_values` ratio (below 1 is faster), 192 cells, **110 faster and none slower beyond 5%**: 32 0.909-0.911, 64 0.873-0.879, 96 0.877, 128 0.789 (1T) to 0.911, 256 0.856 (1T) and 0.999-1.004 threaded, 512 0.922 (1T) and 0.999-1.000 threaded, 1024x64 0.975-0.986, 4096x256 0.978-0.999. The 36 `qr_factor` cells the arm flagged are the binary-layout effect: a control with the kernel compiled but not called moves them by the same amount (`on/off` 1.01-1.20 while `off/pin` is 0.87-1.06). | EXPERIMENT | tlinalg-rs#46 | The predeclared target (`svd_values/64` approaching 0.275 ms) is **not** met (0.47). Three tuning attempts (shift-relative deflation, late-failure shift correction, analytic two-value deflation) each fail to converge on a wide-spectrum block: the missing piece is an admissible shift (`dlasq4`), which is also what 512+ needs. Worklog: `docs/worklogs/2026-10-23-dqds-values-only.md`; evidence `20261023T000823Z-dqds-blast` / `20261023T012045Z-dqds-blast-baseline` / `20261023T043848Z-dqds-off-control`. Tests: `dqds::tests::{matches_the_reference_on_bidiagonals,diagonal_is_its_own_answer,scale_does_not_change_the_answer}`. |
 | blocked bidiagonalisation, the dqds shift, the work copy's leading dimension, and the caller-held worker group (values-only SVD) | `bidiag_in_place` is LAPACK's `dgebrd`/`dlabrd` form: panels of 16 reduced against deferred corrections `X`/`Y`, one trailing update per panel through two plain `matmul` calls, the unblocked loop on the tail below a remaining dimension of 128 (sizes at or below 128 untouched); each panel row is gathered into a contiguous buffer, because a row of `A` has stride `m` and neither the vector products nor the reflector norm take a strided operand. dqds: ping-pong buffers instead of save/restore, release relative to the eigenvalue (`(100 eps)^2 (sigma + a)`), trailing 2x2 in closed form, the shift backed off from the trailing 2x2's eigenvalue by twice its eigenvector-weighted second-order coupling, a rejection at the last factor corrected to `tau + delta`: 3.5 accepted and 0.2 rejected sweeps per value at 512 (from 5.7 and 3.8); values-only takes dqds at every size and budget. The SVD's work copy gets 64 more bytes per column when its stride would be a multiple of 1 KiB (a 4 KiB stride put a row's entries in one L1 set: bidiag 512 8T 18.6 -> 15.3 ms). The worker group is the caller's: the library takes no `spindle::with_lock`, re-exports `faer::spindle` and documents that a caller with a thread budget holds it (the two raw-rayon loops of the block-apply split now go through spindle, as they starved under a held group). Final campaign against the pin under a harness that holds the group (`20261023T125347Z-blocked-dqds-caller-lock-2` vs `20261023T112343Z-caller-lock-baseline`, 384 cells, control 1.0027, **266 faster / 0 slower beyond 5%**): `svd_values` 32-96 0.60-0.63, 128 0.57-0.67, 256 0.60-0.67, 512 0.67 (1T) / 0.61 (8T), 4096x256 0.92 (1T) / 0.57 (8T); `svd_thin` 0.79 at 512 8T, 0.51 at 4096x256 8T. Against MKL in the same run (F order): f64 512 0.85 / 0.91 / 0.97 at 1T / 4T / 8T, 4096x256 0.83 / 0.93 / 1.11, 64 1.17; complex faster at 1T everywhere (0.60-0.99) and at 8T for 512 and both tall shapes. The unblocked loop keeps its row in a contiguous buffer too (`41161e9`: `svd_values` 64 0.296 -> 0.239 ms, 128 1.30 -> 1.01), and a basis of up to eight columns forms its `T` by inner products on the calling thread (`2351819`: the recursive QR's deepest upgrades, 2.0 -> 0.27 ms at 4096x256 8T, `qr_factor` 0.47-0.96 of the pin). Final run `20261023T193712Z-row-buffer-narrow-basis` against same-harness baselines: SVD families 363 faster / 0 slower over 384 cells (control 1.0010), `qr_factor` 174 / 0 over 192 (control 0.9947); `svd_values` 0.51 of the pin at 64, 0.61 at 512 8T, 0.56 at 4096x256 8T; against MKL f64 **64 0.94**, 512 0.85-0.95, 4096x256 0.81 / 0.92 / 1.09. Provenance: published algorithms (Dongarra-Hammarling-Sorensen 1989; Fernando-Parlett 1994; Parlett-Marques LAWN 155) and measured breakdowns only; no vendor or LAPACK code consulted, the vendor timed through its public API. | EXPERIMENT | tlinalg-rs#46 | The predeclared target (`svd_values/64` approaching 0.275 ms) is **met: 0.239 ms against the vendor's 0.255** (`41161e9`). Left: 32 (1.09), the eight-thread tall shape (1.09) and 256 at 8T (1.02), in the QR panel's narrow products and the 256-sized bidiagonalisation that does not thread (`docs/worklogs/2026-10-23-svd-64-and-tall.md`). For the maintainer: adoption of the whole set, and whether `faer::spindle` stays a public re-export. Worklogs: `docs/worklogs/2026-10-23-blocked-bidiag.md`, `docs/worklogs/2026-10-23-dqds-shift.md`; evidence in `faer-benchmark`: `20261023T091426Z-blocked-bidiag`, `20261023T100330Z-blocked-dqds` (old harness, no held group), `20261023T112343Z-caller-lock-baseline`, `20261023T120538Z-blocked-dqds-caller-lock`, `20261023T125347Z-blocked-dqds-caller-lock-2`. Tests: bidiag reconstruction at 200/300/400x300 real and complex, sequential and threaded; `dqds::tests::matches_the_reference_on_a_dense_bidiagonal` (1e-12 per value up to 200); lib 218 passed / 2 ignored. |
 | dqds falls back to the rotation path when it does not finish (#110) | Since the values-only route went to dqds (row above), a bidiagonal whose singular values span more than the range of their squares returned `NoConvergence`. The case is LAPACK's DCHKBD type 16, entries `eps^U(-2, 2)`: the zero-shift sweep's running pivot underflows to zero and every retry is rejected. It hit all four dtypes at n = 127 to 129 in `dchkbd_15_to_19` and `bidiagonal_family`. `compute_bidiag_real_svd` now finishes on the rotation path when `dqds::singular_values` fails. That is LAPACK's structure: `DBDSQR` falls through to its implicit QR when `DLASQ1` reports that dqds did not finish (`INFO = 2`). `dqds::singular_values` writes `d` only on success and never writes `e`, so the fallback starts from the caller's bidiagonal. `DLASQ1`'s scaling (the largest entry to `sqrt(eps / safmin)` before squaring) was tried and does not rescue this matrix, so it is not adopted. Fork-only: the values-only dqds route came from #52 and upstream faer has no dqds, so upstream is not affected. | KEEP | tensor4all/faer#110 | Correctness, no timing claim: only the failing inputs take the second path. Test `dchkbd_type_16_values_only_finishes` (n = 128 at recursion threshold 8 and n = 127 at 128, so both fallback routes; all four `U`/`V` modes, every dtype) runs with the regular tests and fails without the fallback. The opt-in `dchkbd_15_to_19` and `bidiagonal_family` pass in all four dtypes. Record: [worklog](docs/worklogs/2026-10-29-dqds-fallback.md). |
+| dqds takes its work arrays from the caller's `MemStack` (#156) | The values-only SVD allocated on the heap in 0.26.0, though the caller passes a `MemStack` sized by `svd_scratch`; 0.25.0 did not. `dqds::singular_values` built its four work arrays (the squared factors and the ping-pong buffers) with `vec!`, two allocations at 1 x 1 and four above. It also sorted with a stable sort, whose merge buffer goes to the heap once the slice outgrows a stack buffer. The arrays now come from the stack (`dqds::singular_values_scratch`, included in the values-only bidiagonal scratch and so in `svd_scratch`), and the sort is `sort_unstable_by`. Found by tlinalg's `alloc_counts` test, which blocked tlinalg 0.8.0. | KEEP | tensor4all/faer#156 | Test `faer/tests/svd_values_no_alloc.rs`: its own binary with a process-wide counting allocator, values-only `svd` on `Par::Seq` for f64/f32/c64/c32 at 1 x 1, 2 x 2, 17 x 17, 40 x 33, 33 x 40, 128 x 128, 300 x 60 (QR route) and 700 x 700, after a warm-up, counts zero allocations. Before the change it counts `[2, 2, 2, 2]` at 1 x 1 (the issue's numbers); with the work arrays fixed but the stable sort kept, `[1, 0, 1, 0]` at 700 x 700. No arithmetic change. |
 | Full pivoting LU: the pivot search's ties are independent of the SIMD lane count (#102, PR #105) | `reduce_2d`, the SIMD pivot search of `lu::full_pivoting`, took the first lane that held the maximum (`first_true_mask`). Which lane comes first depends on the lane count, so for a matrix with tied maxima (small integers tie a lot) an AVX-512 build and an AVX2 build chose different pivots and eliminated in a different order. On a singular fixture one of them ended on an exact zero pivot, which the solve turned into `NaN`. The tie is now broken as the scalar search breaks it: the smallest column, then the smallest row. Upstream has the same code: faer 0.24.4 `src/linalg/lu/full_pivoting/factor.rs:64` (`simd.first_true_mask(is_best)`). Not reported upstream. | KEEP | tensor4all/faer#102 | Correctness: the pivots are the same on every vector width. Test `pivot_search_breaks_ties_by_column_then_row` (the SIMD and scalar searches agree on tied matrices of several shapes). Upstream report: [`upstream-reports/faer-001-full-piv-lu-lane-tie`](upstream-reports/faer-001-full-piv-lu-lane-tie/). |
 | Partial SVD: a retained Ritz pair carries its residual estimate through the restart (#102, PR #108) | In the restart of the partial SVD (`operator::svd`), a retained pair that moves to a lower slot moved its value and its vector but not its residual estimate. The coupling row of the next iteration is built from the residual array (`H[(i, k)] = residual[i]`), so the moved pair coupled with another pair's estimate. The locked branch already moved it. The retained branch now moves it too. Upstream has the same code: faer 0.24.4 `src/operator/svd/mod.rs:273-275`, whose `Group::Retain` branch moves `H[(mi, mi)]` and the vector, but not `residual`. The partial self-adjoint eigensolver (`operator::self_adjoint_eigen`) is not affected, because it rebuilds the coupling row from `H Q`. Not reported upstream. | KEEP | tensor4all/faer#102 | Correctness. PR #108 adds no dedicated test; its record is the PR. Upstream report: [`upstream-reports/faer-002-partial-svd-residual`](upstream-reports/faer-002-partial-svd-residual/). |
 | Krylov solvers: final residual and BiCGSTAB breakdown (#133, #134; upstream #328, #327) | When `bicgstab` or `conjugate_gradient` stopped at `max_iters`, `NoConvergence` carried the initial residual `|b|`: `bicgstab` kept each new norm in a local, and `conjugate_gradient` shadowed `abs_residual` inside the loop. Both now report the residual of the returned `x`. `bicgstab`'s rank test on `r̃ᵀ v` took its tolerance from its own largest pivot, so a zero first pivot passed and the step divided by zero, turning `x` into NaN until `max_iters`. The tolerance is now `eps * k * max_j |r̃_j| * max_j |v_j|`, the bound on one entry of `r̃ᵀ v`, so a wide block of right-hand sides is not scaled by its width. The comparison is strict, so `v = 0` also counts. At a breakdown (rank 0), or when the half step `s` or the new residual is not finite, the iteration stops. It keeps the last iterate whose residual is finite and reports that residual as `NoConvergence`. A finite residual does not bound `x` itself: an operator scaled near the range limit can still overflow `x`. `InitialGuessStatus::Zero` is now written out as zeros on entry in both solvers, so every exit, including `max_iters = 0`, returns zeros or an iterate, never the caller's buffer. A breakdown-specific error variant is the feature in #135 (upstream #334) and is not part of this change. Written from the upstream reports (Codeberg sarah-quinones/faer#327, #328); no code was taken from the open upstream PRs #330 and #331. | KEEP | tensor4all/faer#133, tensor4all/faer#134 | Correctness only, no timing claim: the rank test adds `2k` column norms of length `n` per iteration, and a zero guess adds one fill. Nine focused tests: `codeberg_328` for both solvers, `codeberg_327`, a first-step breakdown with a NaN-filled zero guess, `A r̃ = 0`, a non-finite second half step, a zero guess with `max_iters = 0` in both solvers, and an `f32` identity with 64 right-hand sides. Before the fix, `codeberg_327`, both `codeberg_328` tests and the first-step breakdown failed. All 23 `matrix_free` tests pass. Independent review: gpt-6.1-sol. |

faer/src/linalg/svd/dqds.rs

@@ -41,13 +41,24 @@
 
 use crate::internal_prelude::*;
 use crate::linalg::svd::SvdError;
-use alloc::vec;
 
 
 /// How many sweeps one entry may take before the block gives up.
 const SWEEPS_PER_ENTRY: usize = 200;
 
-/// Singular values of the bidiagonal with diagonal `d` and off-diagonal `e`, in place.
+/// The workspace of [`singular_values`] for an `n x n` bidiagonal: the squared factors and the
+/// second buffer of the ping-pong, two arrays of `n` and two of `n - 1`.
+pub(crate) fn singular_values_scratch<T: RealField>(n: usize) -> StackReq {
+	StackReq::all_of(&[
+		StackReq::new::<T>(n),
+		StackReq::new::<T>(n.saturating_sub(1)),
+		StackReq::new::<T>(n),
+		StackReq::new::<T>(n.saturating_sub(1)),
+	])
+}
+
+/// Singular values of the bidiagonal with diagonal `d` and off-diagonal `e`, in place. The work
+/// arrays come from `stack` ([`singular_values_scratch`]); nothing is allocated on the heap.
 ///
 /// On success `d` holds the singular values in **descending** order. `e` is never written, and on
 /// an error `d` is unchanged too, so a caller can finish on another path from the same bidiagonal
@@ -57,7 +68,7 @@ const SWEEPS_PER_ENTRY: usize = 200;
 /// `d` must not contain a zero: the transform divides by the factors, and a rank-deficient
 /// bidiagonal has a cheaper answer than this loop. The caller checks that and keeps its own
 /// rotation-based path for it.
-pub(crate) fn singular_values<T: RealField>(d: &mut [T], e: &mut [T]) -> Result<(), SvdError> {
+pub(crate) fn singular_values<T: RealField>(d: &mut [T], e: &mut [T], stack: &mut MemStack) -> Result<(), SvdError> {
     let n = d.len();
     core::assert!(e.len() + 1 >= n, "off-diagonal shorter than the diagonal");
     if n == 0 {
@@ -84,8 +95,8 @@ pub(crate) fn singular_values<T: RealField>(d: &mut [T], e: &mut [T]) -> Result<
     }
 
     let inv_scale = scale.recip();
-    let mut a = vec![zero::<T>(); n];
-    let mut b = vec![zero::<T>(); n.saturating_sub(1)];
+    let (mut a, stack) = stack.make_with::<T>(n, |_| zero());
+    let (mut b, stack) = stack.make_with::<T>(n.saturating_sub(1), |_| zero());
     for i in 0..n {
         let v = d[i].clone() * inv_scale.clone();
         a[i] = v.clone() * v.clone();
@@ -97,8 +108,8 @@ pub(crate) fn singular_values<T: RealField>(d: &mut [T], e: &mut [T]) -> Result<
 
     // The second buffer of the ping-pong: a sweep writes the factors it produces beside the ones it
     // reads, so a rejected sweep costs nothing to undo and an accepted one is a swap of roles.
-    let mut a2 = vec![zero::<T>(); n];
-    let mut b2 = vec![zero::<T>(); n.saturating_sub(1)];
+    let (mut a2, stack) = stack.make_with::<T>(n, |_| zero());
+    let (mut b2, _) = stack.make_with::<T>(n.saturating_sub(1), |_| zero());
 
     // Split at off-diagonals that are negligible against their neighbours before any shift exists.
     let tol2 = {
@@ -124,8 +135,10 @@ pub(crate) fn singular_values<T: RealField>(d: &mut [T], e: &mut [T]) -> Result<
     for i in 0..n {
         d[i] = a[i].clone().fmax(zero::<T>()).sqrt() * scale.clone();
     }
-    // The public contract is descending order; the iteration produces none.
-    d.sort_by(|x, y| y.partial_cmp(x).unwrap_or(core::cmp::Ordering::Equal));
+    // The public contract is descending order; the iteration produces none. Unstable: the stable
+    // sort allocates its merge buffer on the heap once the slice outgrows a stack buffer (700 `f64`
+    // values do).
+    d.sort_unstable_by(|x, y| y.partial_cmp(x).unwrap_or(core::cmp::Ordering::Equal));
     Ok(())
 }
 
@@ -360,6 +373,12 @@ fn two_by_two_shifted<T: RealField>(d1: T, f2: T, d2: T) -> (T, T) {
 mod tests {
     use super::*;
     use crate::Mat;
+    use dyn_stack::MemBuffer;
+
+    /// `singular_values` with its scratch from a fresh buffer
+    fn values<T: RealField>(d: &mut [T], e: &mut [T]) -> Result<(), SvdError> {
+        singular_values(d, e, MemStack::new(&mut MemBuffer::new(singular_values_scratch::<T>(d.len()))))
+    }
 
     /// The kernel against `thin_svd` of the same bidiagonal, over spectra that make the iteration
     /// work: well separated, clustered, repeated, badly scaled, and with negative diagonals.
@@ -413,7 +432,7 @@ mod tests {
 
                 let mut dd = d.clone();
                 let mut ee = e.clone();
-                singular_values(&mut dd, &mut ee).unwrap();
+                values(&mut dd, &mut ee).unwrap();
                 for i in 0..n {
                     let tol = 1e-10 * want[0].max(1.0);
                     core::assert!(
@@ -451,7 +470,7 @@ mod tests {
             bidiag::bidiag_in_place(bid.as_mut(), hl.as_mut(), hr.as_mut(), Par::Seq, MemStack::new(&mut mem), Default::default());
             let mut d: Vec<f64> = (0..n).map(|i| bid[(i, i)]).collect();
             let mut e: Vec<f64> = (0..n - 1).map(|i| bid[(i, i + 1)]).collect();
-            singular_values(&mut d, &mut e).unwrap();
+            values(&mut d, &mut e).unwrap();
             let want: Vec<f64> = a.as_ref().thin_svd(Par::Seq).unwrap().S().column_vector().iter().copied().collect();
             for i in 0..n {
                 core::assert!((d[i] - want[i]).abs() <= 1e-12 * want[i], "n={n} i={i}: {} vs {}", d[i], want[i]);
@@ -464,7 +483,7 @@ mod tests {
     fn diagonal_is_its_own_answer() {
         let mut d = vec![3.0f64, 1.0, 2.0];
         let mut e: Vec<f64> = vec![0.0; 2];
-        singular_values(&mut d, &mut e).unwrap();
+        values(&mut d, &mut e).unwrap();
         core::assert_eq!(d, vec![3.0, 2.0, 1.0]);
     }
 
@@ -475,11 +494,11 @@ mod tests {
         let e = vec![0.3f64, 0.2, 0.1];
         let mut d1 = d.clone();
         let mut e1 = e.clone();
-        singular_values(&mut d1, &mut e1).unwrap();
+        values(&mut d1, &mut e1).unwrap();
 
         let mut d2: Vec<f64> = d.iter().map(|x| x * 1e150).collect();
         let mut e2: Vec<f64> = e.iter().map(|x| x * 1e150).collect();
-        singular_values(&mut d2, &mut e2).unwrap();
+        values(&mut d2, &mut e2).unwrap();
         for i in 0..d1.len() {
             let got = d2[i] * 1e-150;
             core::assert!((got - d1[i]).abs() <= 1e-10 * d1[0], "{got} vs {}", d1[i]);

faer/src/linalg/svd/mod.rs

@@ -166,10 +166,12 @@ fn bidiag_real_svd_scratch<T: RealField>(
 	par: ParShape,
 	params: SvdParams,
 ) -> StackReq {
+	// the values-only route runs dqds first, and falls through to the paths below when it fails
+	let dqds = if compute_u || compute_v { StackReq::EMPTY } else { dqds::singular_values_scratch::<T>(n) };
 	if n < params.recursion_threshold {
-		StackReq::EMPTY
+		dqds
 	} else {
-		StackReq::all_of(&[
+		StackReq::any_of(&[dqds, StackReq::all_of(&[
 			temp_mat_scratch::<T>(2, if compute_u { 0 } else { n + 1 }),
 			bidiag_svd::divide_and_conquer_scratch::<T>(
 				n,
@@ -178,7 +180,7 @@ fn bidiag_real_svd_scratch<T: RealField>(
 				compute_v,
 				par,
 			),
-		])
+		])])
 	}
 }
 fn compute_bidiag_cplx_svd<T: ComplexField>(
@@ -341,7 +343,7 @@ fn compute_bidiag_real_svd<T: RealField>(
 		// units at all (#110). `dqds::singular_values` writes `diag` only on success and never
 		// writes `subdiag`, so the rotation path starts from the caller's bidiagonal.
 		if !diag.rb().iter().any(|x| *x == zero())
-			&& dqds::singular_values(diag.rb_mut().as_slice_mut(), subdiag.rb_mut().as_slice_mut()).is_ok()
+			&& dqds::singular_values(diag.rb_mut().as_slice_mut(), subdiag.rb_mut().as_slice_mut(), stack).is_ok()
 		{
 			return Ok(());
 		}