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

faer/linalg/svd/
mod.rs
