faer/linalg/cholesky/llt/
reconstruct.rs

1use crate::assert;
2use crate::internal_prelude::*;
3use linalg::matmul::triangular::BlockStructure;
4
5pub fn reconstruct_scratch<T: ComplexField>(dim: usize, par: Par) -> StackReq {
6	_ = (dim, par);
7	StackReq::EMPTY
8}
9
10#[track_caller]
11pub fn reconstruct<T: ComplexField>(out: MatMut<'_, T>, L: MatRef<'_, T>, par: Par, stack: &mut MemStack) {
12	let mut out = out;
13	_ = stack;
14
15	let n = out.nrows();
16	assert!(all(out.nrows() == n, out.ncols() == n, L.nrows() == n, L.ncols() == n,));
17
18	linalg::matmul::triangular::matmul(
19		out.rb_mut(),
20		BlockStructure::TriangularLower,
21		Accum::Replace,
22		L,
23		BlockStructure::TriangularLower,
24		L.adjoint(),
25		BlockStructure::TriangularUpper,
26		one(),
27		par,
28	);
29}
30
31#[cfg(test)]
32mod tests {
33	use super::*;
34	use crate::assert;
35	use crate::stats::prelude::*;
36	use crate::utils::approx::*;
37	use dyn_stack::MemBuffer;
38	use linalg::cholesky::llt::*;
39
40	#[test]
41	fn test_reconstruct() {
42		let rng = &mut StdRng::seed_from_u64(0);
43		let n = 50;
44
45		let A = CwiseMatDistribution {
46			nrows: n,
47			ncols: n,
48			dist: ComplexDistribution::new(StandardNormal, StandardNormal),
49		}
50		.rand::<Mat<c64>>(rng);
51
52		let A = &A * A.adjoint();
53		let mut L = A.to_owned();
54
55		factor::cholesky_in_place(
56			L.as_mut(),
57			Default::default(),
58			Par::Seq,
59			MemStack::new(&mut { MemBuffer::new(factor::cholesky_in_place_scratch::<c64>(n, Par::Seq, default())) }),
60			default(),
61		)
62		.unwrap();
63
64		let approx_eq = CwiseMat(ApproxEq::eps() * (n as f64));
65
66		let mut A_rec = Mat::zeros(n, n);
67		reconstruct::reconstruct(
68			A_rec.as_mut(),
69			L.as_ref(),
70			Par::Seq,
71			MemStack::new(&mut MemBuffer::new(reconstruct::reconstruct_scratch::<c64>(n, Par::Seq))),
72		);
73
74		for j in 0..n {
75			for i in 0..j {
76				A_rec[(i, j)] = A_rec[(j, i)].conj();
77			}
78		}
79
80		assert!(A_rec ~ A);
81	}
82}
faer/linalg/cholesky/llt/reconstruct.rs

faer/linalg/cholesky/llt/
reconstruct.rs
