faer/linalg/qr/col_pivoting/
inverse.rs

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

faer/linalg/qr/col_pivoting/
inverse.rs
