lair 0.8.0

Linear algebra in Rust
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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

use std::cmp;
use std::ops::{Div, MulAssign};

use ndarray::{s, Array1, ArrayBase, DataMut, Ix2};

use crate::{lapack, Scalar};

/// Computes the QR factorization of a matrix.
pub fn geqrf<A, S>(a: &mut ArrayBase<S, Ix2>) -> Array1<A>
where
    A: Scalar + Div<<A as Scalar>::Real, Output = A> + MulAssign<<A as Scalar>::Real>,
    S: DataMut<Elem = A>,
{
    let min_dim = cmp::min(a.nrows(), a.ncols());
    let mut tau = Array1::<A>::zeros(min_dim);
    for i in 0..min_dim {
        let bottom = a.nrows();
        let (beta, _, t) = lapack::larfg(a[(i, i)], a.column_mut(i).slice_mut(s![i + 1..bottom]));
        tau[i] = t;
        if i < a.ncols() {
            a[(i, i)] = A::one();
            let v = a.column(i).slice(s![i..]).to_owned();
            lapack::larf::left(&v, t.conj(), &mut a.slice_mut(s![i.., i + 1..]));
            a[(i, i)] = beta.into();
        }
    }
    tau
}

#[cfg(test)]
mod tests {
    use approx::assert_abs_diff_eq;
    use ndarray::{arr1, arr2};
    use num_complex::Complex64;

    #[test]
    fn square_real_smallest() {
        let a = arr2(&[[2_f64]]);
        let mut qr = a.clone();
        assert_eq!(qr.shape(), &[1, 1]);
        let tau = super::geqrf(&mut qr);
        assert_eq!(tau.shape(), &[1]);
        assert_abs_diff_eq!(qr, arr2(&[[2_f64]]), epsilon = 1e-8);
        assert_abs_diff_eq!(tau, arr1(&[0.]), epsilon = 1e-8);
    }

    #[test]
    fn square_real() {
        let a = arr2(&[
            [1_f64, 2_f64, 4_f64],
            [0_f64, 0_f64, 5_f64],
            [0_f64, 3_f64, 6_f64],
        ]);
        let mut qr = a.clone();
        assert_eq!(qr.shape(), &[3, 3]);
        let tau = super::geqrf(&mut qr);
        assert_eq!(tau.shape(), &[3]);
        assert_abs_diff_eq!(
            qr,
            arr2(&[
                [1_f64, 2_f64, 4_f64],
                [0_f64, -3_f64, -6_f64],
                [0_f64, 1_f64, -5_f64],
            ]),
            epsilon = 1e-8
        );
        assert_abs_diff_eq!(tau, arr1(&[0., 1., 0.]), epsilon = 1e-8);
    }

    #[test]
    fn square_complex() {
        let a = arr2(&[
            [Complex64::new(1., 1.), Complex64::new(2., -1.)],
            [Complex64::new(3., 1.), Complex64::new(4., -1.)],
        ]);
        let mut qr = a.clone();
        assert_eq!(qr.shape(), &[2, 2]);
        let tau = super::geqrf(&mut qr);
        assert_eq!(tau.shape(), &[2]);
        assert_abs_diff_eq!(qr[(0, 0)].re, -3.46410162, epsilon = 1e-6);
        assert_abs_diff_eq!(qr[(0, 0)].im, 0., epsilon = 1e-6);
        assert_abs_diff_eq!(qr[(0, 1)].re, -3.46410162, epsilon = 1e-6);
        assert_abs_diff_eq!(qr[(0, 1)].im, 2.88675135, epsilon = 1e-6);
        assert_abs_diff_eq!(qr[(1, 0)].re, 0.68769902, epsilon = 1e-6);
        assert_abs_diff_eq!(qr[(1, 0)].im, 0.0699583, epsilon = 1e-6);
        assert_abs_diff_eq!(qr[(1, 1)].re, 1.29099445, epsilon = 1e-6);
        assert_abs_diff_eq!(qr[(1, 1)].im, 0., epsilon = 1e-6);
        assert_abs_diff_eq!(tau[0].re, 1.28867513, epsilon = 1e-6);
        assert_abs_diff_eq!(tau[0].im, 0.28867513, epsilon = 1e-6);
        assert_abs_diff_eq!(tau[1].re, 1.02290316, epsilon = 1e-6);
        assert_abs_diff_eq!(tau[1].im, -0.99973769, epsilon = 1e-6);
    }

    #[test]
    fn tall() {
        let a = arr2(&[
            [1_f64, 2_f64, 3_f64],
            [2_f64, 2_f64, 1_f64],
            [3_f64, 1_f64, 2_f64],
            [2_f64, 3_f64, 3_f64],
        ]);
        let mut qr = a.clone();
        let tau = super::geqrf(&mut qr);
        assert_abs_diff_eq!(
            qr,
            arr2(&[
                [-4.24264069, -3.53553391, -4.00693843],
                [0.38148714, 2.34520788, 2.06094026],
                [0.57223071, 0.88223561, -1.64224532],
                [0.38148714, -0.36153262, -0.34951992]
            ]),
            epsilon = 1e-6
        );
        assert_abs_diff_eq!(
            tau,
            arr1(&[1.23570226, 1.04764396, 1.7822704]),
            epsilon = 1e-6
        );
    }

    #[test]
    fn wide() {
        let a = arr2(&[
            [1_f32, 2_f32, 3_f32, 1_f32],
            [2_f32, 2_f32, 1_f32, 3_f32],
            [3_f32, 1_f32, 2_f32, 2_f32],
        ]);
        let mut qr = a.clone();
        let tau = super::geqrf(&mut qr);
        assert_abs_diff_eq!(
            qr,
            arr2(&[
                [-3.741_657_5, -2.405_351_2, -2.939_873_7, -3.474_396_2],
                [0.42179344, -1.792_842_9, -1.6334791, -0.916_341_9],
                [0.632_690_2, -0.9237749, -1.639_783_1, 1.043_498_4]
            ]),
            epsilon = 1e-6
        );
        assert_abs_diff_eq!(tau, arr1(&[1.267_261_3, 1.079_121_1, 0.]), epsilon = 1e-6);
    }
}