diff --git a/PATCHSET.md b/PATCHSET.md index 35c9ad66..36c4e3bc 100644 --- a/PATCHSET.md +++ b/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. | diff --git a/faer/src/linalg/svd/dqds.rs b/faer/src/linalg/svd/dqds.rs index adbc0b48..1e9bb139 100644 --- a/faer/src/linalg/svd/dqds.rs +++ b/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(n: usize) -> StackReq { + StackReq::all_of(&[ + StackReq::new::(n), + StackReq::new::(n.saturating_sub(1)), + StackReq::new::(n), + StackReq::new::(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(d: &mut [T], e: &mut [T]) -> Result<(), SvdError> { +pub(crate) fn singular_values(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(d: &mut [T], e: &mut [T]) -> Result< } let inv_scale = scale.recip(); - let mut a = vec![zero::(); n]; - let mut b = vec![zero::(); n.saturating_sub(1)]; + let (mut a, stack) = stack.make_with::(n, |_| zero()); + let (mut b, stack) = stack.make_with::(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(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::(); n]; - let mut b2 = vec![zero::(); n.saturating_sub(1)]; + let (mut a2, stack) = stack.make_with::(n, |_| zero()); + let (mut b2, _) = stack.make_with::(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(d: &mut [T], e: &mut [T]) -> Result< for i in 0..n { d[i] = a[i].clone().fmax(zero::()).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(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(d: &mut [T], e: &mut [T]) -> Result<(), SvdError> { + singular_values(d, e, MemStack::new(&mut MemBuffer::new(singular_values_scratch::(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 = (0..n).map(|i| bid[(i, i)]).collect(); let mut e: Vec = (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 = 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 = 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 = d.iter().map(|x| x * 1e150).collect(); let mut e2: Vec = 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]); diff --git a/faer/src/linalg/svd/mod.rs b/faer/src/linalg/svd/mod.rs index 90db0d91..2caf5766 100644 --- a/faer/src/linalg/svd/mod.rs +++ b/faer/src/linalg/svd/mod.rs @@ -166,10 +166,12 @@ fn bidiag_real_svd_scratch( 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::(n) }; if n < params.recursion_threshold { - StackReq::EMPTY + dqds } else { - StackReq::all_of(&[ + StackReq::any_of(&[dqds, StackReq::all_of(&[ temp_mat_scratch::(2, if compute_u { 0 } else { n + 1 }), bidiag_svd::divide_and_conquer_scratch::( n, @@ -178,7 +180,7 @@ fn bidiag_real_svd_scratch( compute_v, par, ), - ]) + ])]) } } fn compute_bidiag_cplx_svd( @@ -341,7 +343,7 @@ fn compute_bidiag_real_svd( // 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(()); } diff --git a/faer/tests/svd_values_no_alloc.rs b/faer/tests/svd_values_no_alloc.rs new file mode 100644 index 00000000..72cfde2f --- /dev/null +++ b/faer/tests/svd_values_no_alloc.rs @@ -0,0 +1,87 @@ +//! #156: the values-only SVD allocates nothing on the heap once the caller's `MemStack` is sized by +//! `svd_scratch`. The dqds route (0.26.0) built its four work arrays with `vec!` (two allocations at +//! 1 x 1, four above), and sorted with a stable sort, whose merge buffer goes to the heap once the +//! slice outgrows a stack buffer (one allocation for `f64` and `c64` at 700 x 700). The allocator counts every +//! allocation of the process while the call runs; the matrices and the scratch are built outside the +//! counted interval. +#![cfg(feature = "std")] +use faer::dyn_stack::{MemBuffer, MemStack}; +use faer::linalg::svd::{ComputeSvdVectors, svd, svd_scratch}; +use faer::diag::Diag; +use faer::traits::ComplexField; +use faer::{Auto, Mat, Par, ParShape, c32, c64}; +use std::alloc::{GlobalAlloc, Layout, System}; +use std::sync::atomic::{AtomicBool, AtomicUsize, Ordering}; + +struct Counting; +static COUNTING: AtomicBool = AtomicBool::new(false); +static ALLOCATIONS: AtomicUsize = AtomicUsize::new(0); +unsafe impl GlobalAlloc for Counting { + unsafe fn alloc(&self, layout: Layout) -> *mut u8 { + if COUNTING.load(Ordering::Relaxed) { + ALLOCATIONS.fetch_add(1, Ordering::Relaxed); + } + unsafe { System.alloc(layout) } + } + unsafe fn alloc_zeroed(&self, layout: Layout) -> *mut u8 { + if COUNTING.load(Ordering::Relaxed) { + ALLOCATIONS.fetch_add(1, Ordering::Relaxed); + } + unsafe { System.alloc_zeroed(layout) } + } + unsafe fn realloc(&self, ptr: *mut u8, layout: Layout, new_size: usize) -> *mut u8 { + if COUNTING.load(Ordering::Relaxed) { + ALLOCATIONS.fetch_add(1, Ordering::Relaxed); + } + unsafe { System.realloc(ptr, layout, new_size) } + } + unsafe fn dealloc(&self, ptr: *mut u8, layout: Layout) { + unsafe { System.dealloc(ptr, layout) } + } +} +#[global_allocator] +static GLOBAL: Counting = Counting; + +/// Allocations of one values-only `svd` call on an `m x n` matrix, after a warm-up call. +fn values_only_allocations(a: &Mat) -> usize { + let (m, n) = a.shape(); + let params = >::auto().into(); + let mut mem = MemBuffer::new(svd_scratch::(m, n, ComputeSvdVectors::No, ComputeSvdVectors::No, ParShape::Seq, params)); + let mut s = Diag::::zeros(Ord::min(m, n)); + let mut run = || svd(a.as_ref(), s.as_mut(), None, None, Par::Seq, MemStack::new(&mut mem), params).unwrap(); + run(); // warm-up: first-use allocations (lazy statics, feature detection) are not the call's + ALLOCATIONS.store(0, Ordering::SeqCst); + COUNTING.store(true, Ordering::SeqCst); + run(); + COUNTING.store(false, Ordering::SeqCst); + ALLOCATIONS.load(Ordering::SeqCst) +} + +fn sample(m: usize, n: usize, seed: u64) -> Mat { + let mut x = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + Mat::from_fn(m, n, |_, _| { + x = x.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + ((x >> 11) as f64 / (1u64 << 53) as f64) - 0.5 + }) +} + +// one test function: the counter is process wide, so nothing else may run while it counts +#[test] +fn values_only_svd_allocates_nothing() { + // 1 x 1 (tlinalg's case), small, wide, square past the recursion threshold, tall enough to take + // the QR route, and large enough for the stable sort's buffer to leave the stack + for (m, n) in [(1, 1), (2, 2), (17, 17), (40, 33), (33, 40), (128, 128), (300, 60), (700, 700)] { + let re = sample(m, n, 156); + let im = sample(m, n, 157); + let f32s = Mat::::from_fn(m, n, |i, j| re[(i, j)] as f32); + let c64s = Mat::::from_fn(m, n, |i, j| c64::new(re[(i, j)], im[(i, j)])); + let c32s = Mat::::from_fn(m, n, |i, j| c32::new(re[(i, j)] as f32, im[(i, j)] as f32)); + let counts = [ + values_only_allocations(&re), + values_only_allocations(&f32s), + values_only_allocations(&c64s), + values_only_allocations(&c32s), + ]; + assert_eq!(counts, [0; 4], "{m} x {n}: allocations (f64, f32, c64, c32)"); + } +}