faer/linalg/qr/no_pivoting/
reconstruct.rs

1use crate::assert;
2use crate::internal_prelude::*;
3
4pub fn reconstruct_scratch<T: ComplexField>(nrows: usize, ncols: usize, blocksize: usize, par: Par) -> StackReq {
5	_ = par;
6	linalg::householder::apply_block_householder_sequence_on_the_left_in_place_scratch::<T>(nrows, blocksize, ncols)
7}
8
9#[track_caller]
10pub fn reconstruct<T: ComplexField>(
11	out: MatMut<'_, T>,
12	Q_basis: MatRef<'_, T>,
13	Q_coeff: MatRef<'_, T>,
14	R: MatRef<'_, T>,
15	par: Par,
16	stack: &mut MemStack,
17) {
18	let m = Q_basis.nrows();
19	let n = R.ncols();
20	let size = Ord::min(m, n);
21	assert!(all(
22		out.nrows() == m,
23		out.ncols() == n,
24		Q_basis.nrows() == m,
25		Q_basis.ncols() == size,
26		Q_coeff.ncols() == size,
27		R.nrows() == size,
28		R.ncols() == n,
29	));
30
31	let mut out = out;
32	out.fill(zero());
33	out.rb_mut().get_mut(..size, ..n).copy_from_triangular_upper(R);
34
35	linalg::householder::apply_block_householder_sequence_on_the_left_in_place_with_conj(Q_basis, Q_coeff, Conj::No, out.rb_mut(), par, stack);
36}
37
38#[cfg(test)]
39mod tests {
40	use super::*;
41	use crate::assert;
42	use crate::stats::prelude::*;
43	use crate::utils::approx::*;
44	use dyn_stack::MemBuffer;
45	use linalg::qr::no_pivoting::*;
46
47	#[test]
48	fn test_reconstruct() {
49		let rng = &mut StdRng::seed_from_u64(0);
50		for (m, n) in [(100, 50), (50, 100)] {
51			let A = CwiseMatDistribution {
52				nrows: m,
53				ncols: n,
54				dist: ComplexDistribution::new(StandardNormal, StandardNormal),
55			}
56			.rand::<Mat<c64>>(rng);
57			let size = Ord::min(m, n);
58
59			let mut QR = A.to_owned();
60			let mut Q_coeff = Mat::zeros(4, size);
61
62			factor::qr_in_place(
63				QR.as_mut(),
64				Q_coeff.as_mut(),
65				Par::Seq,
66				MemStack::new(&mut { MemBuffer::new(factor::qr_in_place_scratch::<c64>(m, n, 4, Par::Seq, default())) }),
67				default(),
68			);
69
70			let approx_eq = CwiseMat(ApproxEq::eps() * (n as f64));
71
72			let mut A_rec = Mat::zeros(m, n);
73			reconstruct::reconstruct(
74				A_rec.as_mut(),
75				QR.get(.., ..size),
76				Q_coeff.as_ref(),
77				QR.get(..size, ..),
78				Par::Seq,
79				MemStack::new(&mut MemBuffer::new(reconstruct::reconstruct_scratch::<c64>(m, n, 4, Par::Seq))),
80			);
81
82			assert!(A_rec ~ A);
83		}
84	}
85}
faer/linalg/qr/no_pivoting/reconstruct.rs

faer/linalg/qr/no_pivoting/
reconstruct.rs
