faer/linalg/cholesky/bunch_kaufman/
mod.rs

1//! the bunch kaufman decomposition of a self-adjoint matrix $a$ is such that:
2//! $$P A P^\top = LBL^H$$
3//! where $P$ is a permutation matrix, $B$ is a block diagonal matrix, with $1\times 1$ or $2 \times
4//! 2 $ diagonal blocks, and $L$ is a unit lower triangular matrix
5#![allow(missing_docs)]
6
7pub mod factor;
8pub mod solve;
9
10pub mod inverse;
11pub mod reconstruct;
12
13#[cfg(test)]
14mod tests {
15	use super::factor::PivotingStrategy;
16	use super::*;
17	use crate::internal_prelude::*;
18	use crate::stats::prelude::*;
19	use crate::{assert, c64};
20	use dyn_stack::MemBuffer;
21	use factor::LbltParams;
22	use std::vec;
23
24	#[test]
25	fn test_real() {
26		let rng = &mut StdRng::seed_from_u64(0);
27
28		for n in [3, 6, 19, 100, 421] {
29			let a = CwiseMatDistribution {
30				nrows: n,
31				ncols: n,
32				dist: StandardNormal,
33			}
34			.rand::<Mat<f64>>(rng);
35
36			let a = &a + a.adjoint();
37			let rhs = CwiseMatDistribution {
38				nrows: n,
39				ncols: 2,
40				dist: StandardNormal,
41			}
42			.rand::<Mat<f64>>(rng);
43
44			let mut ldl = a.clone();
45			let mut subdiag = Diag::<f64>::zeros(n);
46
47			let mut perm = vec![0usize; n];
48			let mut perm_inv = vec![0; n];
49
50			let params = Default::default();
51			let mut mem = MemBuffer::new(factor::cholesky_in_place_scratch::<usize, f64>(n, Par::Seq, params));
52			let (_, perm) = factor::cholesky_in_place(
53				ldl.as_mut(),
54				subdiag.as_mut(),
55				&mut perm,
56				&mut perm_inv,
57				Par::Seq,
58				MemStack::new(&mut mem),
59				params,
60			);
61
62			let mut mem = MemBuffer::new(solve::solve_in_place_scratch::<usize, f64>(n, rhs.ncols(), Par::Seq));
63			let mut x = rhs.clone();
64			solve::solve_in_place_with_conj(
65				ldl.as_ref(),
66				ldl.diagonal(),
67				subdiag.as_ref(),
68				Conj::No,
69				perm.rb(),
70				x.as_mut(),
71				Par::Seq,
72				MemStack::new(&mut mem),
73			);
74
75			let err = &a * &x - &rhs;
76			let mut max = 0.0;
77			zip!(err.as_ref()).for_each(|unzip!(err)| {
78				let err = err.abs();
79				if err > max {
80					max = err
81				}
82			});
83			assert!(max < 1e-9);
84		}
85	}
86
87	#[test]
88	fn test_cplx() {
89		let rng = &mut StdRng::seed_from_u64(0);
90
91		for n in [2, 3, 6, 19, 100, 421] {
92			let distribution = ComplexDistribution::new(StandardNormal, StandardNormal);
93			let a = CwiseMatDistribution {
94				nrows: n,
95				ncols: n,
96				dist: distribution,
97			}
98			.rand::<Mat<c64>>(rng);
99
100			let A = &a + a.adjoint();
101			let rhs = CwiseMatDistribution {
102				nrows: n,
103				ncols: 2,
104				dist: distribution,
105			}
106			.rand::<Mat<c64>>(rng);
107
108			for pivoting in [
109				PivotingStrategy::Partial,
110				PivotingStrategy::PartialDiag,
111				PivotingStrategy::Rook,
112				PivotingStrategy::RookDiag,
113				PivotingStrategy::Full,
114			] {
115				let mut ldl = A.clone();
116				let mut subdiag = Diag::<c64>::zeros(n);
117
118				let mut perm = vec![0usize; n];
119				let mut perm_inv = vec![0; n];
120
121				let params = LbltParams {
122					pivoting,
123					blocksize: 4,
124					..auto!(c64)
125				};
126				let mut mem = MemBuffer::new(factor::cholesky_in_place_scratch::<usize, c64>(n, Par::Seq, params.into()));
127				let (_, perm) = factor::cholesky_in_place(
128					ldl.as_mut(),
129					subdiag.as_mut(),
130					&mut perm,
131					&mut perm_inv,
132					Par::Seq,
133					MemStack::new(&mut mem),
134					params.into(),
135				);
136
137				let mut x = rhs.clone();
138				let mut mem = MemBuffer::new(solve::solve_in_place_scratch::<usize, c64>(n, rhs.ncols(), Par::Seq));
139				solve::solve_in_place_with_conj(
140					ldl.as_ref(),
141					ldl.diagonal(),
142					subdiag.as_ref(),
143					Conj::Yes,
144					perm.rb(),
145					x.as_mut(),
146					Par::Seq,
147					MemStack::new(&mut mem),
148				);
149
150				let err = A.conjugate() * &x - &rhs;
151				let mut max = 0.0;
152				zip!(err.as_ref()).for_each(|unzip!(err)| {
153					let err = abs(err);
154					if err > max {
155						max = err
156					}
157				});
158				for i in 0..n {
159					assert!(ldl[(i, i)].im == 0.0);
160				}
161				assert!(max < 1e-9);
162			}
163		}
164	}
165}
faer/linalg/cholesky/bunch_kaufman/mod.rs

faer/linalg/cholesky/bunch_kaufman/
mod.rs
