faer/linalg/qr/col_pivoting/
reconstruct.rs

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

faer/linalg/qr/col_pivoting/
reconstruct.rs
