t4a-faer 0.28.1

Dense and sparse linear algebra for HPC applications with caller owned thread teams.
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

#![allow(dead_code, unused_imports)]
pub use dyn_stack::{MemBuffer, MemStack};
pub use faer::traits::ext::*;
pub use faer::traits::math_utils::*;
pub use faer::traits::{ComplexField, RealField};
pub use faer::{Col, Mat, MatRef, Par, ParShape, c32, c64};

pub const SIZES: [usize; 6] = [1, 2, 3, 8, 33, 100];
pub fn sample(i: usize, j: usize) -> f64 {
	// Deterministic pseudo-random entries, independent of the optional rand
	// feature.
	let mut x =
		(i as u64 + 7).wrapping_mul(0x9E3779B97F4A7C15) ^ (j as u64 + 19);
	x ^= x >> 30;
	x = x.wrapping_mul(0xBF58476D1CE4E5B9);
	x ^= x >> 27;
	(x >> 11) as f64 / (1u64 << 53) as f64 - 0.5
}
pub fn tolerance<T: ComplexField>(n: usize) -> T::Real {
	eps::<T::Real>() * from_f64::<T::Real>(200.0 * n.max(1) as f64)
}
pub fn close<T: ComplexField>(a: MatRef<'_, T>, b: MatRef<'_, T>) {
	let error = (a - b).norm_l2();
	let bound = tolerance::<T>(a.nrows()) * b.norm_l2().fmax(one::<T::Real>());
	assert!(error <= bound, "residual {error:?}, bound {bound:?}");
}
pub fn unitary<T: ComplexField>(q: MatRef<'_, T>) {
	close(
		(&q.adjoint() * q).as_ref(),
		Mat::identity(q.ncols(), q.ncols()).as_ref(),
	);
}
pub fn schur_form<T: ComplexField>(t: MatRef<'_, T>, real: bool) {
	let n = t.nrows();
	for j in 0..n {
		for i in j + if real { 2 } else { 1 }..n {
			assert_eq!(t[(i, j)], zero());
		}
		if real && j + 1 < n && t[(j + 1, j)] != zero() {
			assert_eq!(t[(j, j)], t[(j + 1, j + 1)]);
			assert!(t[(j, j + 1)].real() * t[(j + 1, j)].real() < zero());
			if j + 2 < n {
				assert_eq!(t[(j + 2, j + 1)], zero());
			}
		}
	}
}