faer/linalg/svd/
mod.rs

1//! low level implementation of the svd of a matrix
2//!
3//! the svd of a matrix $A$ of shape $(m, n)$ is a decomposition into three
4//! components $U$, $S$, and $V$, such that:
5//!
6//! - $U$ has shape $(m, m)$ and is a unitary matrix
7//! - $V$ has shape $(n, n)$ and is a unitary matrix
8//! - $S$ has shape $(m, n)$ and is zero everywhere except the main diagonal
9//! - and finally:
10//!
11//! $$A = U S V^H$$
12use crate::assert;
13use crate::internal_prelude::*;
14use bidiag::BidiagParams;
15use linalg::qr::no_pivoting::factor::QrParams;
16/// bidiagonalization
17pub mod bidiag;
18pub(crate) mod bidiag_svd;
19pub(crate) mod dqds;
20/// whether the singular vectors should be computed
21#[derive(Copy, Clone, Debug, PartialEq, Eq)]
22pub enum ComputeSvdVectors {
23	/// do not compute singular vectors
24	No,
25	/// compute the first $\min(\text{nrows}, \text{ncols})$ singular vectors
26	Thin,
27	/// compute singular vectors
28	Full,
29}
30/// svd error
31#[derive(Copy, Clone, Debug, PartialEq, Eq)]
32pub enum SvdError {
33	/// reached max iterations
34	NoConvergence,
35}
36/// svd tuning parameters
37#[derive(Debug, Copy, Clone)]
38pub struct SvdParams {
39	/// bidiagonalization parameters
40	pub bidiag: BidiagParams,
41	/// $QR$ parameters
42	pub qr: QrParams,
43	/// threshold at which the implementation should stop recursing
44	pub recursion_threshold: usize,
45	/// threshold at which parallelism should be disabled
46	pub qr_ratio_threshold: f64,
47	#[doc(hidden)]
48	pub non_exhaustive: NonExhaustive,
49}
50impl<T: ComplexField> Auto<T> for SvdParams {
51	fn auto() -> Self {
52		Self {
53			recursion_threshold: 128,
54			qr_ratio_threshold: 11.0 / 6.0,
55			bidiag: auto!(T),
56			qr: auto!(T),
57			non_exhaustive: NonExhaustive(()),
58		}
59	}
60}
61/// rows to allocate for an `m`-row work copy of the input.
62///
63/// a column stride that is a multiple of 1 KiB maps the entries of one row onto a few L1 sets (a 512-row
64/// `f64` copy has a 4 KiB stride and every entry of a row lands in the same set), and the
65/// bidiagonalisation reads and writes rows at every step. 64 more bytes per column spread a row over
66/// the sets; measured on the bidiagonalisation alone at eight threads: 512 18.6 -> 15.3 ms, 256
67/// 5.3 -> 4.8 ms, with 500 and 520 at 14.6 and 15.8.
68fn work_rows<T>(m: usize) -> usize {
69	let size = core::mem::size_of::<T>();
70	if size > 0 && (m * size) % 1024 == 0 {
71		m + 64usize.div_ceil(size)
72	} else {
73		m
74	}
75}
76fn svd_imp_scratch<T: ComplexField>(
77	m: usize,
78	n: usize,
79	compute_u: ComputeSvdVectors,
80	compute_v: ComputeSvdVectors,
81	bidiag_svd_scratch: fn(
82		n: usize,
83		compute_u: bool,
84		compute_v: bool,
85		par: ParShape,
86		params: SvdParams,
87	) -> StackReq,
88	params: SvdParams,
89	par: ParShape,
90) -> StackReq {
91	assert!(m >= n);
92	let householder_block_size =
93		linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
94	let bid = temp_mat_scratch::<T>(work_rows::<T>(m), n);
95	let householder_left = temp_mat_scratch::<T>(householder_block_size, n);
96	let householder_right = temp_mat_scratch::<T>(householder_block_size, n);
97	let compute_bidiag =
98		bidiag::bidiag_in_place_scratch::<T>(m, n, par, params.bidiag.into());
99	let diag = temp_mat_scratch::<T>(n, 1);
100	let subdiag = diag;
101	let compute_ub = compute_v != ComputeSvdVectors::No;
102	let compute_vb = compute_u != ComputeSvdVectors::No;
103	let u_b = temp_mat_scratch::<T>(if compute_ub { n + 1 } else { 2 }, n + 1);
104	let v_b = temp_mat_scratch::<T>(n, if compute_vb { n } else { 0 });
105	let compute_bidiag_svd =
106		bidiag_svd_scratch(n, compute_ub, compute_vb, par, params);
107	let apply_householder_u = linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(
108		m,
109		householder_block_size,
110		match compute_u {
111			ComputeSvdVectors::No => 0,
112			ComputeSvdVectors::Thin => n,
113			ComputeSvdVectors::Full => m,
114		},
115	);
116	let apply_householder_v = linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(
117		n - 1,
118		householder_block_size,
119		match compute_v {
120			ComputeSvdVectors::No => 0,
121			_ => n,
122		},
123	);
124	StackReq::all_of(&[
125		bid,
126		householder_left,
127		householder_right,
128		StackReq::any_of(&[
129			compute_bidiag,
130			StackReq::all_of(&[
131				diag,
132				subdiag,
133				u_b,
134				v_b,
135				StackReq::any_of(&[
136					compute_bidiag_svd,
137					StackReq::all_of(&[
138						apply_householder_u,
139						apply_householder_v,
140					]),
141				]),
142			]),
143		]),
144	])
145}
146fn bidiag_cplx_svd_scratch<T: ComplexField>(
147	n: usize,
148	compute_u: bool,
149	compute_v: bool,
150	par: ParShape,
151	params: SvdParams,
152) -> StackReq {
153	StackReq::all_of(&[
154		temp_mat_scratch::<T>(n, 1).array(4),
155		temp_mat_scratch::<T::Real>(n + 1, if compute_u { n + 1 } else { 0 }),
156		temp_mat_scratch::<T::Real>(n, if compute_v { n } else { 0 }),
157		bidiag_real_svd_scratch::<T::Real>(
158			n, compute_u, compute_v, par, params,
159		),
160	])
161}
162fn bidiag_real_svd_scratch<T: RealField>(
163	n: usize,
164	compute_u: bool,
165	compute_v: bool,
166	par: ParShape,
167	params: SvdParams,
168) -> StackReq {
169	// the values-only route runs dqds first, and falls through to the paths below when it fails
170	let dqds = if compute_u || compute_v { StackReq::EMPTY } else { dqds::singular_values_scratch::<T>(n) };
171	if n < params.recursion_threshold {
172		dqds
173	} else {
174		StackReq::any_of(&[dqds, StackReq::all_of(&[
175			temp_mat_scratch::<T>(2, if compute_u { 0 } else { n + 1 }),
176			bidiag_svd::divide_and_conquer_scratch::<T>(
177				n,
178				params.recursion_threshold,
179				compute_u,
180				compute_v,
181				par,
182			),
183		])])
184	}
185}
186fn compute_bidiag_cplx_svd<T: ComplexField>(
187	mut diag: ColMut<'_, T, usize, ContiguousFwd>,
188	subdiag: ColMut<'_, T, usize, ContiguousFwd>,
189	mut u: Option<MatMut<'_, T>>,
190	mut v: Option<MatMut<'_, T>>,
191	params: SvdParams,
192	par: Par,
193	stack: &mut MemStack,
194) -> Result<(), SvdError> {
195	let n = diag.nrows();
196	let (mut diag_real, stack) =
197		unsafe { temp_mat_uninit::<T::Real, _, _>(n, 1, stack) };
198	let (mut subdiag_real, stack) =
199		unsafe { temp_mat_uninit::<T::Real, _, _>(n, 1, stack) };
200	let (mut u_real, stack) = unsafe {
201		temp_mat_uninit::<T::Real, _, _>(
202			n + 1,
203			if u.is_some() { n + 1 } else { 0 },
204			stack,
205		)
206	};
207	let (mut v_real, stack) = unsafe {
208		temp_mat_uninit::<T::Real, _, _>(
209			n,
210			if v.is_some() { n } else { 0 },
211			stack,
212		)
213	};
214	let (mut col_mul, stack) =
215		unsafe { temp_mat_uninit::<T, _, _>(n, 1, stack) };
216	let (mut row_mul, stack) =
217		unsafe { temp_mat_uninit::<T, _, _>(n - 1, 1, stack) };
218	let mut u_real = u.rb().map(|_| u_real.as_mat_mut());
219	let mut v_real = v.rb().map(|_| v_real.as_mat_mut());
220	let mut diag_real = diag_real
221		.as_mat_mut()
222		.col_mut(0)
223		.try_as_col_major_mut()
224		.unwrap();
225	let mut subdiag_real = subdiag_real
226		.as_mat_mut()
227		.col_mut(0)
228		.try_as_col_major_mut()
229		.unwrap();
230	let mut col_mul = col_mul.as_mat_mut().col_mut(0);
231	let mut row_mul = row_mul.as_mat_mut().col_mut(0);
232	let normalized = |x: T| {
233		if x == zero() {
234			one()
235		} else {
236			let norm1 = x.real().abs().fmax(x.imag().abs());
237			// a subnormal component scale has an infinite reciprocal; rescale it the
238			// way reference `dlarfg` does instead of dividing by a zero product
239			let ref shift = if norm1 >= min_positive::<T::Real>() {
240				from_f64::<T::Real>(1.0)
241			} else {
242				from_f64::<T::Real>(4503599627370496.0)
243			};
244			let y = x.mul_real(shift).mul_real((&norm1 * shift).recip());
245			y.mul_real(y.abs().recip())
246		}
247	};
248	let magnitude = |x: &T| {
249		// `abs` squares components scaled by `sqrt(min_positive)` or
250		// `sqrt(max_positive)`, which flushes a magnitude below `min_positive`
251		// to zero and overflows above about `max_positive / 2`. Outside
252		// `[min_positive, 1 / min_positive]`, rescale by that exact power of
253		// two first; other entries keep the historical `abs` bit for bit.
254		let sml = min_positive::<T::Real>();
255		let big = sml.recip();
256		let m = x.real().abs().fmax(x.imag().abs());
257		if m < sml && m != zero() {
258			x.mul_real(&big).abs() * sml
259		} else if m > big && m.is_finite() {
260			x.mul_real(&sml).abs() * big
261		} else {
262			x.abs()
263		}
264	};
265	let mut col_normalized = normalized(diag[0].conj());
266	col_mul[0] = col_normalized.copy();
267	diag_real[0] = magnitude(&diag[0]);
268	subdiag_real[n - 1] = zero();
269	for i in 1..n {
270		let row_normalized =
271			normalized((&subdiag[i - 1] * col_normalized).conj());
272		subdiag_real[i - 1] = magnitude(&subdiag[i - 1]);
273		row_mul[i - 1] = row_normalized.conj();
274		col_normalized = normalized((&diag[i] * row_normalized).conj());
275		diag_real[i] = magnitude(&diag[i]);
276		col_mul[i] = col_normalized.copy();
277	}
278	compute_bidiag_real_svd(
279		diag_real.rb_mut(),
280		subdiag_real.rb_mut(),
281		u_real.rb_mut(),
282		v_real.rb_mut(),
283		params,
284		par,
285		stack,
286	)?;
287	for i in 0..n {
288		diag[i] = diag_real[i].to_cplx();
289	}
290	let u_real = u_real.rb();
291	let v_real = v_real.rb();
292	if let (Some(mut u), Some(u_real)) = (u.rb_mut(), u_real) {
293		z!(u.rb_mut().row_mut(0), u_real.row(0))
294			.for_each(|uz!(u, r)| *u = r.to_cplx());
295		z!(u.rb_mut().row_mut(n), u_real.row(n))
296			.for_each(|uz!(u, r)| *u = r.to_cplx());
297		for j in 0..u.ncols() {
298			let mut u = u.rb_mut().col_mut(j).subrows_mut(1, n - 1);
299			let u_real = u_real.rb().col(j).subrows(1, n - 1);
300			z!(u.rb_mut(), u_real, row_mul.rb())
301				.for_each(|uz!(u, re, f)| *u = f.mul_real(re));
302		}
303	}
304	if let (Some(mut v), Some(v_real)) = (v.rb_mut(), v_real) {
305		for j in 0..v.ncols() {
306			let mut v = v.rb_mut().col_mut(j);
307			let v_real = v_real.rb().col(j);
308			z!(v.rb_mut(), v_real, col_mul.rb())
309				.for_each(|uz!(v, re, f)| *v = f.mul_real(re));
310		}
311	}
312	Ok(())
313}
314fn compute_bidiag_real_svd<T: RealField>(
315	mut diag: ColMut<'_, T, usize, ContiguousFwd>,
316	mut subdiag: ColMut<'_, T, usize, ContiguousFwd>,
317	mut u: Option<MatMut<'_, T, usize, usize>>,
318	mut v: Option<MatMut<'_, T, usize, usize>>,
319	params: SvdParams,
320	par: Par,
321	stack: &mut MemStack,
322) -> Result<(), SvdError> {
323	let n = diag.nrows();
324	for i in 0..n {
325		if !(diag[i].is_finite() && subdiag[i].is_finite()) {
326			return Err(SvdError::NoConvergence);
327		}
328	}
329	if u.is_none() && v.is_none() {
330		// Values only: dqds applies no rotations, so a sweep costs one division and one
331		// multiplication per entry instead of a rotation's hypot and its two-sided update, and it
332		// is what LAPACK's values-only route uses at every size (`DBDSDC(compq='N')` ->
333		// `DLASDQ` -> `DBDSQR`'s `ROTATE`-false branch -> `DLASQ1`). It is serial, and still
334		// cheaper than the threaded divide-and-conquer at every size and budget measured
335		// (`docs/worklogs/2026-10-23-dqds-shift.md`: 6.6 ms against 13.3 ms at n = 512 on eight
336		// threads, 25 against 43 at 1024). A rank-deficient bidiagonal keeps the rotation path
337		// below, because the qd transform divides by the factors.
338		//
339		// When dqds does not finish, the rotation path below does, as in LAPACK's `DBDSQR`, which
340		// falls through to its implicit QR when `DLASQ1` reports that dqds failed. A bidiagonal
341		// whose singular values span more than the range of their squares (DCHKBD type 16, entries
342		// `eps^U(-2, 2)`: the zero-shift sweep's running pivot underflows) cannot finish in squared
343		// units at all (#110). `dqds::singular_values` writes `diag` only on success and never
344		// writes `subdiag`, so the rotation path starts from the caller's bidiagonal.
345		if !diag.rb().iter().any(|x| *x == zero())
346			&& dqds::singular_values(diag.rb_mut().as_slice_mut(), subdiag.rb_mut().as_slice_mut(), stack).is_ok()
347		{
348			return Ok(());
349		}
350	}
351	if n < params.recursion_threshold {
352		if let Some(mut u) = u.rb_mut() {
353			u.fill(zero());
354			u.diagonal_mut().fill(one());
355		}
356		if let Some(mut v) = v.rb_mut() {
357			v.fill(zero());
358			v.diagonal_mut().fill(one());
359		}
360		bidiag_svd::qr_algorithm(
361			diag.rb_mut(),
362			subdiag.rb_mut(),
363			u.rb_mut().map(|u| u.submatrix_mut(0, 0, n, n)),
364			v.rb_mut(),
365		)?;
366		return Ok(());
367	} else {
368		let (mut u2, stack) = unsafe {
369			temp_mat_uninit::<T::Real, _, _>(
370				2,
371				if u.is_some() { 0 } else { n + 1 },
372				stack,
373			)
374		};
375		bidiag_svd::divide_and_conquer(
376			diag.as_row_shape_mut(n),
377			subdiag.as_row_shape_mut(n),
378			match u {
379				Some(u) => bidiag_svd::MatU::Full(u),
380				None => bidiag_svd::MatU::TwoRowsStorage(u2.as_mat_mut()),
381			},
382			v.map(|m| m.as_shape_mut(n, n)),
383			par,
384			stack,
385			params.recursion_threshold,
386		)
387	}
388}
389/// Workspace for [`bidiagonal_svd`], including the padded
390/// lower-bidiagonal storage used by the existing divide-and-conquer
391/// implementation.
392///
393/// # Examples
394/// ```
395/// use faer::Par;
396/// use faer::linalg::svd::bidiagonal_svd_scratch;
397/// let _buffer = dyn_stack::MemBuffer::new(bidiagonal_svd_scratch::<f64>(
398/// 	3,
399/// 	faer::ParShape::Seq,
400/// 	Default::default(),
401/// ));
402/// ```
403pub fn bidiagonal_svd_scratch<T: ComplexField>(
404	dim: usize,
405	par: ParShape,
406	params: Spec<SvdParams, T>,
407) -> StackReq {
408	let n = dim;
409	if n == 0 {
410		return StackReq::EMPTY;
411	}
412	StackReq::all_of(&[
413		temp_mat_scratch::<T>(n, 1).array(2),
414		temp_mat_scratch::<T>(n + 1, n + 1),
415		temp_mat_scratch::<T>(n, n),
416		bidiag_cplx_svd_scratch::<T>(n, true, true, par, *params),
417	])
418}
419/// Computes a square bidiagonal SVD `B = U * diag(s) * V.adjoint()`.
420/// `diag` is overwritten with nonnegative, descending singular values (with
421/// zero imaginary part). `offdiag` has length `n.saturating_sub(1)` and is
422/// unchanged. `upper` selects the superdiagonal rather than the subdiagonal.
423/// Both `u` and `v` are `n` by `n`; the output is **V**, not its adjoint.
424/// Real and complex entries are supported via the existing phase reduction.
425///
426/// # Panics
427/// Panics on incompatible dimensions or insufficient scratch.
428///
429/// # Errors
430/// Returns [`SvdError::NoConvergence`] for nonfinite input, iteration
431/// failure, or a result that would not be finite (for example when a complex
432/// entry's modulus overflows). Finite input never yields `Ok` with nonfinite
433/// output. Outputs are unspecified on error; `offdiag` is always unchanged.
434///
435/// # Examples
436/// ```
437/// use dyn_stack::{MemBuffer, MemStack};
438/// use faer::linalg::svd::*;
439/// use faer::{Mat, Par};
440/// let mut d = [3.0f64, -2.0];
441/// let (mut u, mut v) = (Mat::zeros(2, 2), Mat::zeros(2, 2));
442/// let mut mem = MemBuffer::new(bidiagonal_svd_scratch::<f64>(
443/// 	2,
444/// 	faer::ParShape::Seq,
445/// 	Default::default(),
446/// ));
447/// bidiagonal_svd(
448/// 	&mut d,
449/// 	&[0.0],
450/// 	true,
451/// 	u.as_mut(),
452/// 	v.as_mut(),
453/// 	Par::Seq,
454/// 	MemStack::new(&mut mem),
455/// 	Default::default(),
456/// )
457/// .unwrap();
458/// assert_eq!(d, [3.0, 2.0]);
459/// ```
460pub fn bidiagonal_svd<T: ComplexField>(
461	diag: &mut [T],
462	offdiag: &[T],
463	upper: bool,
464	mut u: MatMut<'_, T>,
465	mut v: MatMut<'_, T>,
466	par: Par,
467	stack: &mut MemStack,
468	params: Spec<SvdParams, T>,
469) -> Result<(), SvdError> {
470	let n = diag.len();
471	assert!(offdiag.len() == n.saturating_sub(1));
472	assert!(
473		u.nrows() == n && u.ncols() == n && v.nrows() == n && v.ncols() == n
474	);
475	if !diag.iter().chain(offdiag).all(|x| x.is_finite()) {
476		return Err(SvdError::NoConvergence);
477	}
478	if n == 0 {
479		return Ok(());
480	}
481	// SAFETY: Every entry is written below before the kernel reads the
482	// diagonal.
483	let (mut d, stack) = unsafe { temp_mat_uninit::<T, _, _>(n, 1, stack) };
484	let (mut e, stack) = temp_mat_zeroed::<T, _, _>(n, 1, stack);
485	// SAFETY: The phase-reduction kernel fully writes each factor on success;
486	// its error path returns before either factor is read here.
487	let (mut left, stack) =
488		unsafe { temp_mat_uninit::<T, _, _>(n + 1, n + 1, stack) };
489	// SAFETY: Same full-overwrite contract as the left factor above.
490	let (mut right, stack) = unsafe { temp_mat_uninit::<T, _, _>(n, n, stack) };
491	let mut d = d.as_mat_mut().col_mut(0).try_as_col_major_mut().unwrap();
492	let mut e = e.as_mat_mut().col_mut(0).try_as_col_major_mut().unwrap();
493	for i in 0..n {
494		d[i] = if upper {
495			diag[i].conj()
496		} else {
497			diag[i].copy()
498		};
499		if i + 1 < n {
500			e[i] = if upper {
501				offdiag[i].conj()
502			} else {
503				offdiag[i].copy()
504			};
505		}
506	}
507	compute_bidiag_cplx_svd(
508		d.rb_mut(),
509		e.rb_mut(),
510		Some(left.as_mat_mut()),
511		Some(right.as_mat_mut()),
512		*params,
513		par,
514		stack,
515	)?;
516	// Backstop: finite input never reports success with nonfinite output.
517	if !(d.rb().is_all_finite()
518		&& left.as_mat_ref().submatrix(0, 0, n, n).is_all_finite()
519		&& right.as_mat_ref().is_all_finite())
520	{
521		return Err(SvdError::NoConvergence);
522	}
523	for i in 0..n {
524		diag[i] = d[i].copy();
525	}
526	if upper {
527		u.copy_from(right.as_mat_ref());
528		v.copy_from(left.as_mat_ref().submatrix(0, 0, n, n));
529	} else {
530		u.copy_from(left.as_mat_ref().submatrix(0, 0, n, n));
531		v.copy_from(right.as_mat_ref());
532	}
533	Ok(())
534}
535/// bidiag -> divide conquer svd / qr algo
536fn svd_imp<T: ComplexField>(
537	matrix: MatRef<'_, T>,
538	s: ColMut<'_, T>,
539	u: Option<MatMut<'_, T>>,
540	v: Option<MatMut<'_, T>>,
541	bidiag_svd: fn(
542		diag: ColMut<'_, T, usize, ContiguousFwd>,
543		subdiag: ColMut<'_, T, usize, ContiguousFwd>,
544		u: Option<MatMut<'_, T, usize, usize>>,
545		v: Option<MatMut<'_, T, usize, usize>>,
546		params: SvdParams,
547		par: Par,
548		stack: &mut MemStack,
549	) -> Result<(), SvdError>,
550	mut par: Par,
551	stack: &mut MemStack,
552	params: SvdParams,
553) -> Result<(), SvdError> {
554	assert!(matrix.nrows() >= matrix.ncols());
555	let m = matrix.nrows();
556	let n = matrix.ncols();
557	let bs = linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
558	let (mut bid_store, stack) =
559		unsafe { temp_mat_uninit::<T, _, _>(work_rows::<T>(m), n, stack) };
560	let mut bid = bid_store.as_mat_mut().subrows_mut(0, m);
561	let (mut Hl, stack) = unsafe { temp_mat_uninit::<T, _, _>(bs, n, stack) };
562	let (mut Hr, stack) =
563		unsafe { temp_mat_uninit::<T, _, _>(bs, n - 1, stack) };
564	let mut Hl = Hl.as_mat_mut();
565	let mut Hr = Hr.as_mat_mut();
566	bid.copy_from(matrix);
567	bidiag::bidiag_in_place(
568		bid.rb_mut(),
569		Hl.rb_mut(),
570		Hr.rb_mut(),
571		par.reborrow(),
572		stack,
573		params.bidiag.into(),
574	);
575	let (mut diag, stack) = unsafe { temp_mat_uninit::<T, _, _>(n, 1, stack) };
576	let (mut subdiag, stack) =
577		unsafe { temp_mat_uninit::<T, _, _>(n, 1, stack) };
578	let mut diag = diag.as_mat_mut().col_mut(0).try_as_col_major_mut().unwrap();
579	let mut subdiag = subdiag
580		.as_mat_mut()
581		.col_mut(0)
582		.try_as_col_major_mut()
583		.unwrap();
584	let (mut ub, stack) = unsafe {
585		temp_mat_uninit::<T, _, _>(
586			n + 1,
587			if v.is_some() { n + 1 } else { 0 },
588			stack,
589		)
590	};
591	let (mut vb, stack) = unsafe {
592		temp_mat_uninit::<T, _, _>(n, if u.is_some() { n } else { 0 }, stack)
593	};
594	let mut ub = ub.as_mat_mut();
595	let mut vb = vb.as_mat_mut();
596	for i in 0..n {
597		diag[i] = bid[(i, i)].conj();
598		if i + 1 < n {
599			subdiag[i] = bid[(i, i + 1)].conj();
600		} else {
601			subdiag[i] = zero();
602		}
603	}
604	bidiag_svd(
605		diag.rb_mut(),
606		subdiag.rb_mut(),
607		v.rb().map(|_| ub.rb_mut()),
608		u.rb().map(|_| vb.rb_mut()),
609		params,
610		par.reborrow(),
611		stack,
612	)?;
613	{ s }.copy_from(diag);
614	if let Some(mut u) = u {
615		let ncols = u.ncols();
616		u.rb_mut().submatrix_mut(0, 0, n, n).copy_from(vb.rb());
617		u.rb_mut().submatrix_mut(n, 0, m - n, ncols).fill(zero());
618		u.rb_mut().submatrix_mut(0, n, n, ncols - n).fill(zero());
619		u.rb_mut()
620			.submatrix_mut(n, n, ncols - n, ncols - n)
621			.diagonal_mut()
622			.fill(one());
623		linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(bid.rb(), Hl.rb(), Conj::No, u, par.reborrow(), stack);
624	}
625	if let Some(mut v) = v {
626		v.copy_from(ub.rb().submatrix(0, 0, n, n));
627		for j in 1..n {
628			for i in 0..j {
629				bid[(j, i)] = bid[(i, j)].copy();
630			}
631		}
632		linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(
633			bid.rb().submatrix(1, 0, n - 1, n - 1),
634			Hr.rb(),
635			Conj::Yes,
636			v.subrows_mut(1, n - 1),
637			par,
638			stack,
639		);
640	}
641	Ok(())
642}
643fn compute_squareish_svd<T: ComplexField>(
644	matrix: MatRef<'_, T>,
645	s: ColMut<'_, T>,
646	u: Option<MatMut<'_, T>>,
647	v: Option<MatMut<'_, T>>,
648	par: Par,
649	stack: &mut MemStack,
650	params: SvdParams,
651) -> Result<(), SvdError> {
652	if const { T::IS_REAL } {
653		svd_imp::<T::Real>(
654			unsafe { core::mem::transmute(matrix) },
655			unsafe { core::mem::transmute(s) },
656			unsafe { core::mem::transmute(u) },
657			unsafe { core::mem::transmute(v) },
658			compute_bidiag_real_svd::<T::Real>,
659			par,
660			stack,
661			params,
662		)
663	} else {
664		svd_imp::<T>(
665			matrix,
666			s,
667			u,
668			v,
669			compute_bidiag_cplx_svd::<T>,
670			par,
671			stack,
672			params,
673		)
674	}
675}
676/// computes the layout of the workspace required to compute a matrix's svd
677pub fn svd_scratch<T: ComplexField>(
678	nrows: usize,
679	ncols: usize,
680	compute_u: ComputeSvdVectors,
681	compute_v: ComputeSvdVectors,
682	par: ParShape,
683	params: Spec<SvdParams, T>,
684) -> StackReq {
685	let params = params.config;
686	let mut m = nrows;
687	let mut n = ncols;
688	let mut compute_u = compute_u;
689	let mut compute_v = compute_v;
690	if n > m {
691		core::mem::swap(&mut m, &mut n);
692		core::mem::swap(&mut compute_u, &mut compute_v);
693	}
694	if n == 0 {
695		return StackReq::EMPTY;
696	}
697	let bidiag_svd_scratch = if const { T::IS_REAL } {
698		bidiag_real_svd_scratch::<T::Real>
699	} else {
700		bidiag_cplx_svd_scratch::<T>
701	};
702	if m as f64 / n as f64 <= params.qr_ratio_threshold {
703		svd_imp_scratch::<T>(
704			m,
705			n,
706			compute_u,
707			compute_v,
708			bidiag_svd_scratch,
709			params,
710			par,
711		)
712	} else {
713		let bs =
714			linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
715		StackReq::all_of(&[
716			temp_mat_scratch::<T>(m, n),
717			temp_mat_scratch::<T>(bs, n),
718			StackReq::any_of(&[
719				StackReq::all_of(&[
720					temp_mat_scratch::<T>(n, n),
721					svd_imp_scratch::<T>(n, n, compute_u, compute_v, bidiag_svd_scratch, params, par),
722				]),
723				linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(
724					m,
725					bs,
726					match compute_u {
727						ComputeSvdVectors::No => 0,
728						ComputeSvdVectors::Thin => n,
729						ComputeSvdVectors::Full => m,
730					},
731				),
732			]),
733		])
734	}
735}
736/// computes the svd of $A$, with the singular vectors being omitted, thin or
737/// full
738///
739/// the singular are stored in $S$, and the singular vectors in $U$ and $V$ such
740/// that the singular values are sorted in nonincreasing order
741pub fn svd<T: ComplexField>(
742	A: MatRef<'_, T>,
743	s: DiagMut<'_, T>,
744	u: Option<MatMut<'_, T>>,
745	v: Option<MatMut<'_, T>>,
746	mut par: Par,
747	stack: &mut MemStack,
748	params: Spec<SvdParams, T>,
749) -> Result<(), SvdError> {
750	let params = params.config;
751	let (m, n) = A.shape();
752	let size = Ord::min(m, n);
753	assert!(s.dim() == size);
754	let s = s.column_vector_mut();
755	if let Some(u) = u.rb() {
756		assert!(all(
757			u.nrows() == A.nrows(),
758			any(u.ncols() == A.nrows(), u.ncols() == size),
759		));
760	}
761	if let Some(v) = v.rb() {
762		assert!(all(
763			v.nrows() == A.ncols(),
764			any(v.ncols() == A.ncols(), v.ncols() == size),
765		));
766	}
767	#[cfg(feature = "perf-warn")]
768	match (u.rb(), v.rb()) {
769		(Some(matrix), _) | (_, Some(matrix)) => {
770			if matrix.row_stride().unsigned_abs() != 1
771				&& crate::__perf_warn!(QR_WARN)
772			{
773				if matrix.col_stride().unsigned_abs() == 1 {
774					log::warn!(
775						target : "faer_perf",
776						"SVD prefers column-major singular vector matrices. Found row-major matrix."
777					);
778				} else {
779					log::warn!(
780						target : "faer_perf",
781						"SVD prefers column-major singular vector matrices. Found matrix with generic strides."
782					);
783				}
784			}
785		},
786		_ => {},
787	}
788	let mut u = u;
789	let mut v = v;
790	let mut matrix = A;
791	let do_transpose = n > m;
792	if do_transpose {
793		matrix = matrix.transpose();
794		core::mem::swap(&mut u, &mut v)
795	}
796	let (m, n) = matrix.shape();
797	if n == 0 {
798		if let Some(mut u) = u {
799			u.fill(zero());
800			u.rb_mut().diagonal_mut().fill(one());
801		}
802		return Ok(());
803	}
804	if m as f64 / n as f64 <= params.qr_ratio_threshold {
805		compute_squareish_svd(
806			matrix,
807			s,
808			u.rb_mut(),
809			v.rb_mut(),
810			par,
811			stack,
812			params,
813		)?;
814	} else {
815		let bs =
816			linalg::qr::no_pivoting::factor::recommended_block_size::<T>(m, n);
817		let (mut qr, stack) =
818			unsafe { temp_mat_uninit::<T, _, _>(m, n, stack) };
819		let mut qr = qr.as_mat_mut();
820		let (mut householder, stack) =
821			unsafe { temp_mat_uninit::<T, _, _>(bs, n, stack) };
822		let mut householder = householder.as_mat_mut();
823		{
824			qr.copy_from(matrix.rb());
825			linalg::qr::no_pivoting::factor::qr_in_place(
826				qr.rb_mut(),
827				householder.rb_mut(),
828				par.reborrow(),
829				stack,
830				params.qr.into(),
831			);
832		}
833		{
834			let (mut r, stack) =
835				unsafe { temp_mat_uninit::<T, _, _>(n, n, stack) };
836			let mut r = r.as_mat_mut();
837			z!(r.rb_mut()).for_each_triangular_lower(
838				linalg::zip::Diag::Skip,
839				|uz!(dst)| *dst = zero(),
840			);
841			z!(r.rb_mut(), qr.rb().submatrix(0, 0, n, n))
842				.for_each_triangular_upper(
843					linalg::zip::Diag::Include,
844					|uz!(dst, src)| *dst = src.copy(),
845				);
846			compute_squareish_svd(
847				r.rb(),
848				s,
849				u.rb_mut().map(|u| u.submatrix_mut(0, 0, n, n)),
850				v.rb_mut(),
851				par.reborrow(),
852				stack,
853				params,
854			)?;
855		}
856		if let Some(mut u) = u.rb_mut() {
857			u.rb_mut().subrows_mut(n, m - n).fill(zero());
858			if u.ncols() == m {
859				// the extra columns are `Q * [0; I]`, so their top block has
860				// to start at zero instead of at the caller's old entries
861				u.rb_mut().submatrix_mut(0, n, n, m - n).fill(zero());
862				u.rb_mut()
863					.submatrix_mut(n, n, m - n, m - n)
864					.diagonal_mut()
865					.fill(one());
866			}
867			linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(
868				qr.rb(),
869				householder.rb(),
870				Conj::No,
871				u.rb_mut(),
872				par,
873				stack,
874			);
875		}
876	}
877	if do_transpose {
878		if let Some(u) = u.rb_mut() {
879			z!(u).for_each(|uz!(u)| *u = u.conj())
880		}
881		if let Some(v) = v.rb_mut() {
882			z!(v).for_each(|uz!(v)| *v = v.conj())
883		}
884	}
885	Ok(())
886}
887/// computes the layout of the workspace required to compute a matrix's
888/// pseudoinverse, given the svd
889pub fn pseudoinverse_from_svd_scratch<T: ComplexField>(
890	nrows: usize,
891	ncols: usize,
892	par: ParShape,
893) -> StackReq {
894	_ = par;
895	let size = Ord::min(nrows, ncols);
896	StackReq::all_of(&[
897		temp_mat_scratch::<T>(nrows, size),
898		temp_mat_scratch::<T>(ncols, size),
899	])
900}
901/// computes a self-adjoint matrix's pseudoinverse, given the svd factors $S$,
902/// $U$ and $V$
903pub fn pseudoinverse_from_svd<T: ComplexField>(
904	pinv: MatMut<'_, T>,
905	s: DiagRef<'_, T>,
906	u: MatRef<'_, T>,
907	v: MatRef<'_, T>,
908	par: Par,
909	stack: &mut MemStack,
910) {
911	pseudoinverse_from_svd_with_tolerance(
912		pinv,
913		s,
914		u,
915		v,
916		zero(),
917		eps::<T::Real>()
918			* from_f64::<T::Real>(Ord::max(u.nrows(), v.nrows()) as f64),
919		par,
920		stack,
921	);
922}
923/// computes a self-adjoint matrix's pseudoinverse, given the svd factors $S$,
924/// $U$ and $V$, and tolerance parameters for determining zero singular values
925pub fn pseudoinverse_from_svd_with_tolerance<T: ComplexField>(
926	pinv: MatMut<'_, T>,
927	s: DiagRef<'_, T>,
928	u: MatRef<'_, T>,
929	v: MatRef<'_, T>,
930	abs_tol: T::Real,
931	rel_tol: T::Real,
932	par: Par,
933	stack: &mut MemStack,
934) {
935	let mut pinv = pinv;
936	let m = u.nrows();
937	let n = v.nrows();
938	let size = Ord::min(m, n);
939	assert!(all(
940		u.nrows() == m,
941		v.nrows() == n,
942		u.ncols() >= size,
943		v.ncols() >= size,
944		s.dim() >= size
945	));
946	let s = s.column_vector();
947	let u = u.get(.., ..size);
948	let v = v.get(.., ..size);
949	let smax = s.norm_max();
950	let tol = abs_tol.fmax(rel_tol * smax);
951	let (mut u_trunc, stack) =
952		unsafe { temp_mat_uninit::<T, _, _>(m, size, stack) };
953	let (mut vp_trunc, _) =
954		unsafe { temp_mat_uninit::<T, _, _>(n, size, stack) };
955	let mut u_trunc = u_trunc.as_mat_mut();
956	let mut vp_trunc = vp_trunc.as_mat_mut();
957	let mut len = 0;
958	for j in 0..size {
959		let x = s[j].absmax();
960		if x > tol {
961			let ref p = s[j].real().recip();
962			u_trunc.rb_mut().col_mut(len).copy_from(u.col(j));
963			z!(vp_trunc.rb_mut().col_mut(len), v.col(j))
964				.for_each(|uz!(dst, src)| *dst = src.mul_real(p));
965			len += 1;
966		}
967	}
968	let u_trunc = u_trunc.get(.., ..len);
969	let vp_trunc = vp_trunc.get(.., ..len);
970	linalg::matmul::matmul(
971		pinv.rb_mut(),
972		Accum::Replace,
973		vp_trunc,
974		u_trunc.adjoint(),
975		one(),
976		par,
977	);
978}
979#[cfg(test)]
980mod tests {
981	mod qr_stall;
982	// Instantiate the same cases and checks at each native precision.
983	mod lapack_svd_f64 {
984		type Real = f64;
985		type Scalar = Real;
986		const GUARD: Scalar = -12345.0;
987		const ROTATION_PHASE: Scalar = 1.0;
988		fn diagonal_phase(i: usize) -> Scalar {
989			if i % 2 == 0 { 1.0 } else { -1.0 }
990		}
991		use super::super::{
992			bidiag_real_svd_scratch as bidiagonal_scratch,
993			compute_bidiag_real_svd as compute_bidiagonal,
994		};
995		mod lapack {
996			include!("tests/lapack.rs");
997		}
998		mod svd_views {
999			include!("tests/svd_views.rs");
1000		}
1001	}
1002	mod lapack_svd_f32 {
1003		type Real = f32;
1004		type Scalar = Real;
1005		const GUARD: Scalar = -12345.0;
1006		const ROTATION_PHASE: Scalar = 1.0;
1007		fn diagonal_phase(i: usize) -> Scalar {
1008			if i % 2 == 0 { 1.0 } else { -1.0 }
1009		}
1010		use super::super::{
1011			bidiag_real_svd_scratch as bidiagonal_scratch,
1012			compute_bidiag_real_svd as compute_bidiagonal,
1013		};
1014		mod lapack {
1015			include!("tests/lapack.rs");
1016		}
1017		mod svd_views {
1018			include!("tests/svd_views.rs");
1019		}
1020		#[test]
1021		fn f32_fixtures_use_single_precision_limits() {
1022			let (_, large) = lapack::diagonal_case(2, 2, 6);
1023			assert_eq!(large[0] as f64, f32::MAX.sqrt() as f64);
1024			let (_, small) = lapack::diagonal_case(2, 2, 7);
1025			assert_eq!(small[0] as f64, f32::MIN_POSITIVE.sqrt() as f64);
1026			let (_, geometric) = lapack::diagonal_case(3, 3, 4);
1027			assert_eq!(geometric[2] as f64, f32::EPSILON as f64);
1028			let (_, cluster) = lapack::diagonal_case(4, 4, 9);
1029			assert_eq!(cluster[3] as f64, (1.0 - 24.0 * f32::EPSILON) as f64);
1030		}
1031	}
1032	mod lapack_svd_c64 {
1033		type Real = f64;
1034		type Scalar = crate::c64;
1035		const GUARD: Scalar = Scalar::new(-12345.0, 6789.0);
1036		const ROTATION_PHASE: Scalar = Scalar::new(0.0, 1.0);
1037		fn diagonal_phase(i: usize) -> Scalar {
1038			match i % 4 {
1039				0 => Scalar::new(0.0, 1.0),
1040				1 => Scalar::new(0.6, 0.8),
1041				2 => Scalar::new(-0.8, 0.6),
1042				_ => Scalar::new(-1.0, 0.0),
1043			}
1044		}
1045		use super::super::{
1046			bidiag_cplx_svd_scratch as bidiagonal_scratch,
1047			compute_bidiag_cplx_svd as compute_bidiagonal,
1048		};
1049		mod lapack {
1050			include!("tests/lapack.rs");
1051		}
1052		mod svd_views {
1053			include!("tests/svd_views.rs");
1054		}
1055		mod complex_tests {
1056			include!("tests/complex.rs");
1057		}
1058	}
1059	mod lapack_svd_c32 {
1060		type Real = f32;
1061		type Scalar = crate::c32;
1062		const GUARD: Scalar = Scalar::new(-12345.0, 6789.0);
1063		const ROTATION_PHASE: Scalar = Scalar::new(0.0, 1.0);
1064		fn diagonal_phase(i: usize) -> Scalar {
1065			match i % 4 {
1066				0 => Scalar::new(0.0, 1.0),
1067				1 => Scalar::new(0.6, 0.8),
1068				2 => Scalar::new(-0.8, 0.6),
1069				_ => Scalar::new(-1.0, 0.0),
1070			}
1071		}
1072		use super::super::{
1073			bidiag_cplx_svd_scratch as bidiagonal_scratch,
1074			compute_bidiag_cplx_svd as compute_bidiagonal,
1075		};
1076		mod lapack {
1077			include!("tests/lapack.rs");
1078		}
1079		mod svd_views {
1080			include!("tests/svd_views.rs");
1081		}
1082		mod complex_tests {
1083			include!("tests/complex.rs");
1084		}
1085		#[test]
1086		fn c32_fixtures_use_single_precision_limits() {
1087			let (_, large) = lapack::diagonal_case(2, 2, 6);
1088			assert_eq!(large[0] as f64, f32::MAX.sqrt() as f64);
1089			let (_, small) = lapack::diagonal_case(2, 2, 7);
1090			assert_eq!(small[0] as f64, f32::MIN_POSITIVE.sqrt() as f64);
1091			let (_, geometric) = lapack::diagonal_case(3, 3, 4);
1092			assert_eq!(geometric[2] as f64, f32::EPSILON as f64);
1093			let (_, cluster) = lapack::diagonal_case(4, 4, 9);
1094			assert_eq!(cluster[3] as f64, (1.0 - 24.0 * f32::EPSILON) as f64);
1095		}
1096	}
1097	use super::*;
1098	use crate::assert;
1099	use crate::stats::prelude::*;
1100	use crate::utils::approx::*;
1101	use dyn_stack::MemBuffer;
1102	#[track_caller]
1103	fn test_svd<T: ComplexField>(mat: MatRef<'_, T>) {
1104		let (m, n) = mat.shape();
1105		let params = Spec::new(SvdParams {
1106			recursion_threshold: 8,
1107			qr_ratio_threshold: 1.0,
1108			..auto!(T)
1109		});
1110		let approx_eq = CwiseMat(
1111			ApproxEq::<T::Real>::eps()
1112				* from_f64::<T::Real>(8.0 * Ord::max(m, n) as f64).sqrt(),
1113		);
1114		{
1115			let mut s = Mat::zeros(m, n);
1116			let mut u = Mat::zeros(m, m);
1117			let mut v = Mat::zeros(n, n);
1118			svd(
1119				mat.as_ref(),
1120				s.as_mut().diagonal_mut(),
1121				Some(u.as_mut()),
1122				Some(v.as_mut()),
1123				Par::Seq,
1124				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1125					m,
1126					n,
1127					ComputeSvdVectors::Full,
1128					ComputeSvdVectors::Full,
1129					ParShape::Seq,
1130					params,
1131				))),
1132				params,
1133			)
1134			.unwrap();
1135			let reconstructed = &u * &s * v.adjoint();
1136			assert!(reconstructed ~ mat);
1137		}
1138		let size = Ord::min(m, n);
1139		let mut s = Mat::zeros(size, size);
1140		let mut u = Mat::zeros(m, size);
1141		let mut v = Mat::zeros(n, size);
1142		{
1143			svd(
1144				mat.as_ref(),
1145				s.as_mut().diagonal_mut(),
1146				Some(u.as_mut()),
1147				Some(v.as_mut()),
1148				Par::Seq,
1149				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1150					m,
1151					n,
1152					ComputeSvdVectors::Thin,
1153					ComputeSvdVectors::Thin,
1154					ParShape::Seq,
1155					params,
1156				))),
1157				params,
1158			)
1159			.unwrap();
1160			let reconstructed = &u * &s * v.adjoint();
1161			assert!(reconstructed ~ mat);
1162		}
1163		{
1164			let mut s2 = Mat::zeros(size, size);
1165			let mut u2 = Mat::zeros(m, size);
1166			svd(
1167				mat.as_ref(),
1168				s2.as_mut().diagonal_mut(),
1169				Some(u2.as_mut()),
1170				None,
1171				Par::Seq,
1172				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1173					m,
1174					n,
1175					ComputeSvdVectors::Thin,
1176					ComputeSvdVectors::No,
1177					ParShape::Seq,
1178					params,
1179				))),
1180				params,
1181			)
1182			.unwrap();
1183			assert!(s2 ~ s);
1184			assert!(u2 ~ u);
1185		}
1186		{
1187			let mut s2 = Mat::zeros(size, size);
1188			let mut v2 = Mat::zeros(n, size);
1189			svd(
1190				mat.as_ref(),
1191				s2.as_mut().diagonal_mut(),
1192				None,
1193				Some(v2.as_mut()),
1194				Par::Seq,
1195				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1196					m,
1197					n,
1198					ComputeSvdVectors::No,
1199					ComputeSvdVectors::Thin,
1200					ParShape::Seq,
1201					params,
1202				))),
1203				params,
1204			)
1205			.unwrap();
1206			assert!(s2 ~ s);
1207			assert!(v2 ~ v);
1208		}
1209		{
1210			let mut s2 = Mat::zeros(size, size);
1211			svd(
1212				mat.as_ref(),
1213				s2.as_mut().diagonal_mut(),
1214				None,
1215				None,
1216				Par::Seq,
1217				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1218					m,
1219					n,
1220					ComputeSvdVectors::No,
1221					ComputeSvdVectors::No,
1222					ParShape::Seq,
1223					params,
1224				))),
1225				params,
1226			)
1227			.unwrap();
1228			assert!(s2 ~ s);
1229		}
1230	}
1231	fn check_tiny_split<T: ComplexField>(a: MatRef<'_, T>) {
1232		let (m, n) = a.shape();
1233		let size = Ord::min(m, n);
1234		let compute = |threshold, want_u, want_v| {
1235			let params = Spec::new(SvdParams {
1236				recursion_threshold: threshold,
1237				..auto!(T)
1238			});
1239			let mut s = Col::zeros(size);
1240			let mut u = Mat::zeros(m, if want_u { size } else { 0 });
1241			let mut v = Mat::zeros(n, if want_v { size } else { 0 });
1242			let vectors = |want| {
1243				if want {
1244					ComputeSvdVectors::Thin
1245				} else {
1246					ComputeSvdVectors::No
1247				}
1248			};
1249			svd(
1250				a,
1251				s.as_mut().as_diagonal_mut(),
1252				if want_u { Some(u.as_mut()) } else { None },
1253				if want_v { Some(v.as_mut()) } else { None },
1254				Par::Seq,
1255				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1256					m,
1257					n,
1258					vectors(want_u),
1259					vectors(want_v),
1260					ParShape::Seq,
1261					params,
1262				))),
1263				params,
1264			)
1265			.unwrap();
1266			(s, u, v)
1267		};
1268		let (reference, _, _) = compute(usize::MAX, true, true);
1269		for (want_u, want_v) in
1270			[(true, true), (false, false), (true, false), (false, true)]
1271		{
1272			let (s, u, v) = compute(128, want_u, want_v);
1273			assert!((&s - &reference).norm_max() < from_f64::<T::Real>(1e-13));
1274			if want_u {
1275				assert!(
1276					(u.adjoint() * &u - Mat::<T>::identity(size, size))
1277						.norm_l2() < from_f64::<T::Real>(1e-13)
1278				);
1279			}
1280			if want_v {
1281				assert!(
1282					(v.adjoint() * &v - Mat::<T>::identity(size, size))
1283						.norm_l2() < from_f64::<T::Real>(1e-13)
1284				);
1285			}
1286			if want_u && want_v {
1287				let back = &u * s.as_diagonal() * v.adjoint();
1288				assert!(
1289					(&back - a).norm_l2() / a.norm_l2()
1290						< from_f64::<T::Real>(1e-13)
1291				);
1292			}
1293		}
1294	}
1295	#[test]
1296	fn test_v_only_rank_one_exact_scratch() {
1297		fn check<T: ComplexField>(entry: T, tol: f64) {
1298			let n = 8;
1299			let a = Mat::<T>::from_fn(n, n, |i, j| {
1300				if i == 0 && j == 0 {
1301					entry.copy()
1302				} else {
1303					zero()
1304				}
1305			});
1306			let mut params = Spec::<SvdParams, T>::default();
1307			params.recursion_threshold = 4;
1308			let mut s = Col::<T>::zeros(n);
1309			let mut v = Mat::<T>::zeros(n, n);
1310			let mut mem = MemBuffer::new(svd_scratch::<T>(
1311				n,
1312				n,
1313				ComputeSvdVectors::No,
1314				ComputeSvdVectors::Thin,
1315				ParShape::Seq,
1316				params,
1317			));
1318			svd(
1319				a.as_ref(),
1320				s.as_mut().as_diagonal_mut(),
1321				None,
1322				Some(v.as_mut()),
1323				Par::Seq,
1324				MemStack::new(&mut mem),
1325				params,
1326			)
1327			.unwrap();
1328			assert!(s.is_all_finite() && v.is_all_finite());
1329			let expected =
1330				Col::<T>::from_fn(n, |i| if i == 0 { one() } else { zero() });
1331			let tol = from_f64::<T::Real>(tol);
1332			assert!((&s - &expected).norm_max() <= tol.copy());
1333			assert!(
1334				(v.adjoint() * &v - Mat::<T>::identity(n, n)).norm_max() <= tol
1335			);
1336		}
1337		check::<f32>(1.0, 1e-5);
1338		check::<c32>(c32::new(0.0, 1.0), 1e-5);
1339		check::<f64>(1.0, 1e-13);
1340		check::<c64>(c64::new(0.0, 1.0), 1e-13);
1341	}
1342	#[test]
1343	fn test_tiny_diagonal_at_split() {
1344		for n in [128, 256] {
1345			for tiny in [0.0, 1e-15, 1e-13, 1e-10, 1e-6] {
1346				let a = Mat::<f64>::from_fn(n, n, |i, j| {
1347					if (i == j && i < 64) || (j == i + 1 && i < 64) {
1348						1.0
1349					} else if i == 64 && j == 64 {
1350						tiny
1351					} else {
1352						0.0
1353					}
1354				});
1355				check_tiny_split(a.as_ref());
1356				let complex =
1357					Mat::<c64>::from_fn(n, n, |i, j| c64::new(0.0, a[(i, j)]));
1358				check_tiny_split(complex.as_ref());
1359			}
1360		}
1361	}
1362	/// A zero block at a divide-and-conquer split must still hand its parent merge the
1363	/// identity boundary entry that carries the nonzero part of the spectrum, whatever
1364	/// singular vectors the caller asked for (#33).
1365	#[test]
1366	fn test_zero_block_at_split() {
1367		fn check<T: ComplexField>(
1368			n: usize,
1369			threshold: usize,
1370			want: (bool, bool),
1371			tol: f64,
1372		) {
1373			let (want_u, want_v) = want;
1374			let a = Mat::<T>::from_fn(n, n, |i, j| {
1375				if i == n / 2 && j == n / 2 {
1376					one()
1377				} else {
1378					zero()
1379				}
1380			});
1381			let params = Spec::new(SvdParams {
1382				recursion_threshold: threshold,
1383				..auto!(T)
1384			});
1385			let vectors = |want| {
1386				if want {
1387					ComputeSvdVectors::Full
1388				} else {
1389					ComputeSvdVectors::No
1390				}
1391			};
1392			let mut s = Col::<T>::zeros(n);
1393			let mut u = Mat::<T>::zeros(n, if want_u { n } else { 0 });
1394			let mut v = Mat::<T>::zeros(n, if want_v { n } else { 0 });
1395			svd(
1396				a.as_ref(),
1397				s.as_mut().as_diagonal_mut(),
1398				if want_u { Some(u.as_mut()) } else { None },
1399				if want_v { Some(v.as_mut()) } else { None },
1400				Par::Seq,
1401				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1402					n,
1403					n,
1404					vectors(want_u),
1405					vectors(want_v),
1406					ParShape::Seq,
1407					params,
1408				))),
1409				params,
1410			)
1411			.unwrap();
1412			// `norm_max` skips NaN entries, so finiteness is stated separately
1413			assert!(s.is_all_finite());
1414			if want_u {
1415				assert!(u.is_all_finite());
1416			}
1417			if want_v {
1418				assert!(v.is_all_finite());
1419			}
1420			let expected =
1421				Col::<T>::from_fn(n, |i| if i == 0 { one() } else { zero() });
1422			let tol = from_f64::<T::Real>(tol);
1423			assert!((&s - &expected).norm_max() <= tol.copy());
1424			let identity = Mat::<T>::identity(n, n);
1425			if want_u && want_v {
1426				let back = &u * s.as_diagonal() * v.adjoint();
1427				assert!((&back - &a).norm_max() <= tol.copy());
1428			}
1429			if want_u {
1430				assert!(
1431					(u.adjoint() * &u - &identity).norm_max() <= tol.copy()
1432				);
1433			}
1434			if want_v {
1435				assert!(
1436					(v.adjoint() * &v - &identity).norm_max() <= tol.copy()
1437				);
1438			}
1439		}
1440		for n in [8usize, 16, 32] {
1441			for threshold in [4, 8] {
1442				for want in
1443					[(true, true), (false, false), (true, false), (false, true)]
1444				{
1445					check::<f32>(n, threshold, want, 1e-5);
1446					check::<c32>(n, threshold, want, 1e-5);
1447					check::<f64>(n, threshold, want, 1e-13);
1448					check::<c64>(n, threshold, want, 1e-13);
1449				}
1450			}
1451		}
1452		// the reported `singular_values()` call runs at the default threshold
1453		let threshold = <SvdParams as Auto<f64>>::auto().recursion_threshold;
1454		check::<f64>(256, threshold, (false, false), 1e-13);
1455	}
1456	/// An extreme but finite scale must not overflow or underflow the reduction that
1457	/// bidiagonalisation and the SVD share (#36): the spectrum and the reconstruction
1458	/// stay finite and accurate.
1459	#[test]
1460	fn test_extreme_scale() {
1461		fn check<T: ComplexField>(s: f64, tol: f64) {
1462			let scale = from_f64::<T>(s);
1463			let a = Mat::<T>::from_fn(2, 2, |i, j| match (i, j) {
1464				(0, 0) => scale.copy(),
1465				(0, 1) => scale.mul_real(from_f64::<T::Real>(0.5)),
1466				(1, 1) => scale.copy(),
1467				_ => zero(),
1468			});
1469			let params = Spec::new(SvdParams {
1470				recursion_threshold: 4,
1471				..auto!(T)
1472			});
1473			let mut sigma = Mat::<T>::zeros(2, 2);
1474			let mut u = Mat::<T>::zeros(2, 2);
1475			let mut v = Mat::<T>::zeros(2, 2);
1476			svd(
1477				a.as_ref(),
1478				sigma.as_mut().diagonal_mut(),
1479				Some(u.as_mut()),
1480				Some(v.as_mut()),
1481				Par::Seq,
1482				MemStack::new(&mut MemBuffer::new(svd_scratch::<T>(
1483					2,
1484					2,
1485					ComputeSvdVectors::Thin,
1486					ComputeSvdVectors::Thin,
1487					ParShape::Seq,
1488					params,
1489				))),
1490				params,
1491			)
1492			.unwrap();
1493			let rec = &u * &sigma * v.adjoint();
1494			let err = (&rec - &a).norm_max();
1495			assert!(rec.is_all_finite());
1496			assert!(err <= a.norm_max() * from_f64::<T::Real>(tol));
1497		}
1498		for s in [1.0f64, 1e160, 1e-300] {
1499			check::<f64>(s, 1e-14);
1500			check::<c64>(s, 1e-14);
1501		}
1502		for s in [1.0f64, 1e30, 1e-35] {
1503			check::<f32>(s, 1e-6);
1504			check::<c32>(s, 1e-6);
1505		}
1506	}
1507	/// A rank-one matrix's bidiagonalisation cascades down to subnormal entries, where the
1508	/// row norm's reciprocal (and, for complex input, the phase normalisation's) overflows.
1509	#[test]
1510	fn test_rank_one_subnormal_cascade() {
1511		for n in [16usize, 32, 150] {
1512			test_svd(Mat::<f32>::full(n, n, 1.0).as_ref());
1513			test_svd(Mat::<c32>::full(n, n, c32::ONE).as_ref());
1514		}
1515		for n in [32usize, 150] {
1516			test_svd(Mat::<f64>::full(n, n, 1.0).as_ref());
1517			test_svd(Mat::<c64>::full(n, n, c64::ONE).as_ref());
1518		}
1519	}
1520	#[test]
1521	fn test_zero_child_at_split() {
1522		for n in [128, 256] {
1523			for pos in [0, n / 2, n - 1] {
1524				let a = Mat::<f64>::from_fn(n, n, |i, j| {
1525					if i == pos && j == pos { 1.0 } else { 0.0 }
1526				});
1527				check_tiny_split(a.as_ref());
1528				let complex =
1529					Mat::<c64>::from_fn(n, n, |i, j| c64::new(0.0, a[(i, j)]));
1530				check_tiny_split(complex.as_ref());
1531			}
1532		}
1533	}
1534	#[test]
1535	fn test_real() {
1536		let rng = &mut StdRng::seed_from_u64(1);
1537		for (m, n) in [
1538			(3, 2),
1539			(2, 2),
1540			(4, 4),
1541			(15, 10),
1542			(10, 10),
1543			(15, 15),
1544			(50, 50),
1545			(100, 100),
1546			(150, 150),
1547			(150, 20),
1548			(20, 150),
1549		] {
1550			let mat = CwiseMatDistribution {
1551				nrows: m,
1552				ncols: n,
1553				dist: StandardNormal,
1554			}
1555			.rand::<Mat<f64>>(rng);
1556			test_svd(mat.as_ref());
1557		}
1558	}
1559	#[test]
1560	fn test_cplx() {
1561		let rng = &mut StdRng::seed_from_u64(1);
1562		for (m, n) in [
1563			(1, 1),
1564			(2, 2),
1565			(3, 2),
1566			(2, 2),
1567			(3, 3),
1568			(4, 4),
1569			(15, 10),
1570			(10, 10),
1571			(15, 15),
1572			(16, 16),
1573			(17, 17),
1574			(18, 18),
1575			(19, 19),
1576			(20, 20),
1577			(30, 30),
1578			(50, 50),
1579			(100, 100),
1580			(150, 150),
1581			(150, 20),
1582			(20, 150),
1583		] {
1584			let mat = CwiseMatDistribution {
1585				nrows: m,
1586				ncols: n,
1587				dist: ComplexDistribution::new(StandardNormal, StandardNormal),
1588			}
1589			.rand::<Mat<c64>>(rng);
1590			test_svd(mat.as_ref());
1591		}
1592	}
1593	#[test]
1594	fn test_special() {
1595		for (m, n) in [
1596			(3, 2),
1597			(2, 2),
1598			(4, 4),
1599			(15, 10),
1600			(10, 10),
1601			(15, 15),
1602			(50, 50),
1603			(100, 100),
1604			(150, 150),
1605			(150, 20),
1606			(20, 150),
1607		] {
1608			test_svd(Mat::<f64>::zeros(m, n).as_ref());
1609			test_svd(Mat::<c64>::zeros(m, n).as_ref());
1610			test_svd(Mat::<f64>::full(m, n, 1.0).as_ref());
1611			test_svd(Mat::<c64>::full(m, n, c64::ONE).as_ref());
1612			test_svd(Mat::<f64>::identity(m, n).as_ref());
1613			test_svd(Mat::<c64>::identity(m, n).as_ref());
1614		}
1615	}
1616	#[test]
1617	fn test_zink() {
1618		let diag = [
1619			-9.931833701529301,
1620			-10.920807536026027,
1621			-52.33647796311243,
1622			2.3685025127736967,
1623			2.421701994236093,
1624			-0.5051763005624579,
1625			-0.04808263896606017,
1626			-0.003875251886338955,
1627			-0.0006413264967716465,
1628			-0.003381944152463707,
1629			2.981152313236375e-5,
1630			5.4290648208388795e-6,
1631			-6.329275972084404e-7,
1632			-6.879142344209158e-7,
1633			-5.265228263479126e-9,
1634			-2.941999902335516e-9,
1635			-1.3060984997930294e-10,
1636			7.07516117218088e-12,
1637			1.8657003929029376e-12,
1638			-6.216080089659131e-14,
1639		];
1640		let subdiag = [
1641			-57.8029649868477,
1642			17.67263066467847,
1643			8.884153814270894,
1644			-9.01998231080713,
1645			-1.028638150814966,
1646			0.22247719217200435,
1647			0.016389886745811315,
1648			-0.004090989452162578,
1649			0.00036818904090536926,
1650			-0.0031394146217732367,
1651			-7.571300829706796e-6,
1652			3.0045718718618155e-6,
1653			2.1329796886727743e-6,
1654			9.259701025627789e-8,
1655			2.2291214755992877e-9,
1656			-2.3017207713252894e-9,
1657			6.807967994979358e-11,
1658			2.1677299575405587e-12,
1659			-3.07282771050034e-13,
1660			0.0,
1661		];
1662		let n = diag.len();
1663		let params = SvdParams {
1664			recursion_threshold: 8,
1665			qr_ratio_threshold: 1.0,
1666			..auto!(f64)
1667		};
1668		let mut d = ColRef::from_slice(&diag).to_owned();
1669		let mut s = ColRef::from_slice(&subdiag).to_owned();
1670		compute_bidiag_real_svd(
1671			d.as_mut().try_as_col_major_mut().unwrap(),
1672			s.as_mut().try_as_col_major_mut().unwrap(),
1673			None,
1674			None,
1675			params,
1676			Par::Seq,
1677			MemStack::new(&mut MemBuffer::new(bidiag_real_svd_scratch::<f64>(
1678				n,
1679				false,
1680				false,
1681				ParShape::Seq,
1682				params,
1683			))),
1684		)
1685		.unwrap();
1686		assert!(d[n - 1] != 0.0);
1687	}
1688	#[test]
1689	fn test_pinv() {
1690		let rng = &mut StdRng::seed_from_u64(0);
1691		let m = 6;
1692		let n = 36;
1693		let mat = CwiseMatDistribution {
1694			nrows: m,
1695			ncols: n,
1696			dist: StandardNormal,
1697		}
1698		.rand::<Mat<f64>>(rng);
1699		let pinv = mat.svd(Par::Seq).unwrap().pseudoinverse(Par::Seq);
1700		let err = &mat * &pinv - Mat::<f64>::identity(m, m);
1701		assert!(err.norm_max() < 1e-10);
1702	}
1703	#[test]
1704	fn test_pinv2() {
1705		let rng = &mut StdRng::seed_from_u64(0);
1706		let m = 36;
1707		let n = 6;
1708		let mat = CwiseMatDistribution {
1709			nrows: m,
1710			ncols: n,
1711			dist: StandardNormal,
1712		}
1713		.rand::<Mat<f64>>(rng);
1714		let pinv = mat.svd(Par::Seq).unwrap().pseudoinverse(Par::Seq);
1715		let err = &pinv * &mat - Mat::<f64>::identity(n, n);
1716		assert!(err.norm_max() < 1e-10);
1717	}
1718	/// The phase reduction's magnitude and phase helpers outside the normal
1719	/// range: for an exact power of two `s`, the singular values of `s B` are
1720	/// `s` times those of `B`. Entries below `min_positive` were flushed to
1721	/// zero or made `inf`/NaN, entries above `1 / min_positive` overflowed.
1722	fn cplx_phase_reduction_is_scale_exact<T: ComplexField>(
1723		make: fn(f64, f64) -> T,
1724	) {
1725		let n = 3;
1726		let base_d = [make(1.8, -0.4), make(-0.3, 0.7), make(0.5, 0.5)];
1727		let base_e = [make(0.2, 0.6), make(-0.8, 0.1)];
1728		let params: SvdParams = auto!(T);
1729		let run = |s: &T::Real| {
1730			let mut d = Col::<T>::from_fn(n, |i| base_d[i].mul_real(s));
1731			let mut e = Col::<T>::from_fn(n, |i| {
1732				if i + 1 < n {
1733					base_e[i].mul_real(s)
1734				} else {
1735					zero()
1736				}
1737			});
1738			let mut u = Mat::<T>::zeros(n + 1, n + 1);
1739			let mut v = Mat::<T>::zeros(n, n);
1740			compute_bidiag_cplx_svd(
1741				d.as_mut().try_as_col_major_mut().unwrap(),
1742				e.as_mut().try_as_col_major_mut().unwrap(),
1743				Some(u.as_mut()),
1744				Some(v.as_mut()),
1745				params,
1746				Par::Seq,
1747				MemStack::new(&mut MemBuffer::new(
1748					bidiag_cplx_svd_scratch::<T>(
1749						n,
1750						true,
1751						true,
1752						ParShape::Seq,
1753						params,
1754					),
1755				)),
1756			)
1757			.unwrap();
1758			assert!(u.is_all_finite() && v.is_all_finite());
1759			d
1760		};
1761		let reference = run(&one());
1762		let sml = min_positive::<T::Real>();
1763		for s in [sml.recip(), &sml * from_f64::<T::Real>(0.0625)] {
1764			let d = run(&s);
1765			for i in 0..n {
1766				let scaled_back = d[i].real() / &s;
1767				assert!(
1768					(scaled_back - reference[i].real()).abs()
1769						<= from_f64::<T::Real>(1e3)
1770							* eps::<T::Real>() * reference[0].real()
1771				);
1772			}
1773		}
1774	}
1775	#[test]
1776	fn cplx_phase_reduction_outside_normal_range() {
1777		cplx_phase_reduction_is_scale_exact::<c64>(c64::new);
1778		cplx_phase_reduction_is_scale_exact::<c32>(|r, i| {
1779			c32::new(r as f32, i as f32)
1780		});
1781		cplx_phase_reduction_is_scale_exact::<f64>(|r, _| r);
1782		cplx_phase_reduction_is_scale_exact::<f32>(|r, _| r as f32);
1783	}
1784	/// `svd` promises nonnegative singular values in nonincreasing order, for any
1785	/// finite input. This 128 x 128 matrix with three nonzero entries
1786	/// bidiagonalises to a cascade whose spectrum reaches the subnormal range; the
1787	/// divide-and-conquer merge used to classify a nonzero subnormal value with the
1788	/// deflated zeros and return `[0.596..., 1.9e-310, 0, ..., 6.3e-4, ...]` (#82).
1789	#[test]
1790	fn test_sorted_with_subnormal_at_split() {
1791		let n = 128;
1792		// (row, column, real part, imaginary part) of each nonzero entry, as
1793		// reported in #82.
1794		let entries = [
1795			(41, 80, -1.3368711811580607e-64, 3.48452997681476e-64),
1796			(67, 80, -5.871970588424443e-1, -1.065196368760547e-1),
1797			(5, 112, -2.8153454475366503e-4, 5.652970997533713e-4),
1798		];
1799		let mut a = Mat::<c64>::zeros(n, n);
1800		for (row, col, re, im) in entries {
1801			a[(row, col)] = c64::new(re, im);
1802		}
1803		let svd = a.as_ref().svd(Par::Seq).unwrap();
1804		let s = svd.S().column_vector();
1805		let values: std::vec::Vec<f64> =
1806			(0..s.nrows()).map(|i| s[i].re).collect();
1807		for i in 0..n {
1808			assert!(values[i] >= 0.0);
1809			assert!(s[i].im == 0.0);
1810			if i + 1 < n {
1811				assert!(
1812					values[i] >= values[i + 1],
1813					"s[{i}] = {:e} < s[{}] = {:e}",
1814					values[i],
1815					i + 1,
1816					values[i + 1],
1817				);
1818			}
1819		}
1820		// the two nonzero values of the input, as `zgesdd` reports them
1821		assert!((values[0] - 0.596780377486914).abs() <= 1e-14);
1822		assert!((values[1] - 6.315239590698253e-4).abs() <= 1e-17);
1823		let u = svd.U();
1824		let v = svd.V();
1825		let rec = &u * svd.S() * v.adjoint();
1826		assert!((&rec - &a).norm_max() <= 1e-14);
1827		assert!(
1828			(u.adjoint() * &u - Mat::<c64>::identity(n, n)).norm_max() <= 1e-14
1829		);
1830		assert!(
1831			(v.adjoint() * &v - Mat::<c64>::identity(n, n)).norm_max() <= 1e-14
1832		);
1833	}
1834	/// full factors must be orthogonal even when the caller's output buffers
1835	/// hold old entries; the tall/wide QR path used to leave the top-right block
1836	/// of `U` (and, after transposition, of `V`) as it found it
1837	#[test]
1838	fn full_factors_ignore_old_output_entries() {
1839		for (m, n) in [(4, 2), (2, 4)] {
1840			let a = Mat::<f64>::identity(m, n);
1841			let size = Ord::min(m, n);
1842			let mut s = Diag::<f64>::zeros(size);
1843			let mut u = Mat::<f64>::full(m, m, 1.0);
1844			let mut v = Mat::<f64>::full(n, n, 1.0);
1845			svd(
1846				a.as_ref(),
1847				s.as_mut(),
1848				Some(u.as_mut()),
1849				Some(v.as_mut()),
1850				Par::Seq,
1851				MemStack::new(&mut MemBuffer::new(svd_scratch::<f64>(
1852					m,
1853					n,
1854					ComputeSvdVectors::Full,
1855					ComputeSvdVectors::Full,
1856					ParShape::Seq,
1857					default(),
1858				))),
1859				default(),
1860			)
1861			.unwrap();
1862			let orth_u =
1863				(u.transpose() * &u - Mat::<f64>::identity(m, m)).norm_max();
1864			let orth_v =
1865				(v.transpose() * &v - Mat::<f64>::identity(n, n)).norm_max();
1866			assert!(orth_u < 1e-12, "U orthogonality error {orth_u}");
1867			assert!(orth_v < 1e-12, "V orthogonality error {orth_v}");
1868			let rec = u.as_ref().get(.., ..size) * s.as_ref()
1869				* v.as_ref().get(.., ..size).transpose();
1870			assert!((&rec - &a).norm_max() < 1e-12);
1871		}
1872	}
1873}
faer/linalg/svd/mod.rs

faer/linalg/svd/
mod.rs
