-
Notifications
You must be signed in to change notification settings - Fork 65
/
householder.rs
99 lines (90 loc) · 2.51 KB
/
householder.rs
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
use ndarray::*;
use ndarray_linalg::{krylov::*, *};
fn over<A: Scalar + Lapack>(rtol: A::Real) {
const N: usize = 4;
let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
let a: Array2<A> = random_using((N, N * 2), &mut rng);
// Terminate
let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
let a_sub = a.slice(s![.., 0..N]);
let qc: Array2<A> = conjugate(&q);
assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I");
assert_close_l2!(&q.dot(&r), &a_sub, rtol; "Check A = QR");
// Skip
let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip);
let a_sub = a.slice(s![.., 0..N]);
let qc: Array2<A> = conjugate(&q);
assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol);
assert_close_l2!(&q.dot(&r), &a_sub, rtol);
// Full
let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Full);
let qc: Array2<A> = conjugate(&q);
assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol);
assert_close_l2!(&q.dot(&r), &a, rtol);
}
#[test]
fn over_f32() {
over::<f32>(1e-5);
}
#[test]
fn over_f64() {
over::<f64>(1e-9);
}
#[test]
fn over_c32() {
over::<c32>(1e-5);
}
#[test]
fn over_c64() {
over::<c64>(1e-9);
}
fn full<A: Scalar + Lapack>(rtol: A::Real) {
const N: usize = 5;
let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
let a: Array2<A> = random_using((N, N), &mut rng);
let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
let qc: Array2<A> = conjugate(&q);
assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I");
assert_close_l2!(&q.dot(&r), &a, rtol; "Check A = QR");
}
#[test]
fn full_f32() {
full::<f32>(1e-5);
}
#[test]
fn full_f64() {
full::<f64>(1e-9);
}
#[test]
fn full_c32() {
full::<c32>(1e-5);
}
#[test]
fn full_c64() {
full::<c64>(1e-9);
}
fn half<A: Scalar + Lapack>(rtol: A::Real) {
const N: usize = 4;
let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
let a: Array2<A> = random_using((N, N / 2), &mut rng);
let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
let qc: Array2<A> = conjugate(&q);
assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol; "Check Q^H Q = I");
assert_close_l2!(&q.dot(&r), &a, rtol; "Check A = QR");
}
#[test]
fn half_f32() {
half::<f32>(1e-5);
}
#[test]
fn half_f64() {
half::<f64>(1e-9);
}
#[test]
fn half_c32() {
half::<c32>(1e-5);
}
#[test]
fn half_c64() {
half::<c64>(1e-9);
}