faer/linalg/qr/no_pivoting/
inverse.rs

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

faer/linalg/qr/no_pivoting/
inverse.rs
