faer 0.24.0

linear algebra library
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

use crate::assert;
use crate::internal_prelude::*;
use linalg::matmul::triangular::BlockStructure;
pub fn reconstruct_scratch<I: Index, T: ComplexField>(
	nrows: usize,
	ncols: usize,
	par: Par,
) -> StackReq {
	_ = par;
	temp_mat_scratch::<T>(nrows, ncols)
}
#[track_caller]
pub fn reconstruct<I: Index, T: ComplexField>(
	out: MatMut<'_, T>,
	L: MatRef<'_, T>,
	U: MatRef<'_, T>,
	row_perm: PermRef<'_, I>,
	par: Par,
	stack: &mut MemStack,
) {
	let m = L.nrows();
	let n = U.ncols();
	let size = Ord::min(m, n);
	assert!(all(out.nrows() == m, out.ncols() == n, row_perm.len() == m));
	let (mut tmp, _) = unsafe { temp_mat_uninit::<T, _, _>(m, n, stack) };
	let mut tmp = tmp.as_mat_mut();
	let mut out = out;
	linalg::matmul::triangular::matmul(
		tmp.rb_mut().get_mut(..size, ..size),
		BlockStructure::Rectangular,
		Accum::Replace,
		L.get(..size, ..size),
		BlockStructure::UnitTriangularLower,
		U.get(..size, ..size),
		BlockStructure::TriangularUpper,
		one(),
		par,
	);
	if m > n {
		linalg::matmul::triangular::matmul(
			tmp.rb_mut().get_mut(size.., ..size),
			BlockStructure::Rectangular,
			Accum::Replace,
			L.get(size.., ..size),
			BlockStructure::Rectangular,
			U.get(..size, ..size),
			BlockStructure::TriangularUpper,
			one(),
			par,
		);
	}
	if m < n {
		linalg::matmul::triangular::matmul(
			tmp.rb_mut().get_mut(..size, size..),
			BlockStructure::Rectangular,
			Accum::Replace,
			L.get(..size, ..size),
			BlockStructure::UnitTriangularLower,
			U.get(..size, size..),
			BlockStructure::Rectangular,
			one(),
			par,
		);
	}
	crate::perm::permute_rows(out.rb_mut(), tmp.rb(), row_perm.inverse());
}
#[cfg(test)]
mod tests {
	use super::*;
	use crate::assert;
	use crate::stats::prelude::*;
	use crate::utils::approx::*;
	use dyn_stack::MemBuffer;
	use linalg::lu::partial_pivoting::*;
	#[test]
	fn test_reconstruct() {
		let rng = &mut StdRng::seed_from_u64(0);
		for (m, n) in [(100, 50), (50, 100)] {
			let A = CwiseMatDistribution {
				nrows: m,
				ncols: n,
				dist: ComplexDistribution::new(StandardNormal, StandardNormal),
			}
			.rand::<Mat<c64>>(rng);
			let mut LU = A.to_owned();
			let perm_fwd = &mut *vec![0usize; m];
			let perm_bwd = &mut *vec![0usize; m];
			let (_, perm) = factor::lu_in_place(
				LU.as_mut(),
				perm_fwd,
				perm_bwd,
				Par::Seq,
				MemStack::new(&mut {
					MemBuffer::new(factor::lu_in_place_scratch::<usize, c64>(
						m,
						n,
						Par::Seq,
						default(),
					))
				}),
				default(),
			);
			let approx_eq = CwiseMat(ApproxEq::eps() * (n as f64));
			let mut A_rec = Mat::zeros(m, n);
			reconstruct::reconstruct(
				A_rec.as_mut(),
				LU.as_ref(),
				LU.as_ref(),
				perm,
				Par::Seq,
				MemStack::new(&mut MemBuffer::new(
					reconstruct::reconstruct_scratch::<usize, c64>(
						m,
						n,
						Par::Seq,
					),
				)),
			);
			assert!(A_rec ~ A);
		}
	}
}