nalgebra/tests/linalg/ldl.rs at 00ca97b0df453d64b2d06f9eda6c04a7c3404c3b · dimforge/nalgebra · GitHub

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
use na::{Complex, Matrix3};
use num::Zero;

#[test]
#[rustfmt::skip]
fn ldl_simple() {
    let m = Matrix3::new(
        Complex::new(2.0, 0.0), Complex::new(-1.0, 0.5),  Complex::zero(),
        Complex::new(-1.0, -0.5),  Complex::new(2.0, 0.0), Complex::new(-1.0, 0.0),
        Complex::zero(), Complex::new(-1.0, 0.0),  Complex::new(2.0, 0.0));

    let ldl = m.lower_triangle().ldl().unwrap();

    // Rebuild
    let p = ldl.l * ldl.d_matrix() * ldl.l.adjoint();

    assert!(relative_eq!(m, p, epsilon = 3.0e-12));
}

#[test]
#[rustfmt::skip]
fn ldl_partial() {
    let m = Matrix3::new(
        Complex::new(2.0, 0.0), Complex::zero(),  Complex::zero(),
        Complex::zero(),  Complex::zero(), Complex::zero(),
        Complex::zero(), Complex::zero(),  Complex::new(2.0, 0.0));

    let ldl = m.lower_triangle().ldl().unwrap();

    // Rebuild
    let p = ldl.l * ldl.d_matrix() * ldl.l.adjoint();

    assert!(relative_eq!(m, p, epsilon = 3.0e-12));
}

#[test]
#[rustfmt::skip]
fn ldl_cholesky() {
    let m = Matrix3::new(
        Complex::new(2.0, 0.0), Complex::new(-1.0, 0.5),  Complex::zero(),
        Complex::new(-1.0, -0.5),  Complex::new(2.0, 0.0), Complex::new(-1.0, 0.0),
        Complex::zero(), Complex::new(-1.0, 0.0),  Complex::new(2.0, 0.0));

    let chol= m.cholesky().unwrap();
    let ldl = m.ldl().unwrap();

    assert!(relative_eq!(ldl.cholesky_l(), chol.l(), epsilon = 3.0e-16));
}

#[test]
#[should_panic]
#[rustfmt::skip]
fn ldl_non_sym_panic() {
    let m = Matrix3::new(
        2.0, -1.0,  0.0,
        1.0, -2.0,  3.0,
       -2.0,  1.0,  0.3);

    let ldl = m.ldl().unwrap();

    // Rebuild
    let p = ldl.l * ldl.d_matrix() * ldl.l.transpose();

    assert!(relative_eq!(m, p, epsilon = 3.0e-16));
}

#[cfg(feature = "proptest-support")]
mod proptest_tests {
    #[allow(unused_imports)]
    use crate::core::helper::{RandComplex, RandScalar};

    macro_rules! gen_tests(
        ($module: ident, $scalar: expr) => {
            mod $module {
                #[allow(unused_imports)]
                use crate::core::helper::{RandScalar, RandComplex};
                use crate::proptest::*;
                use proptest::{prop_assert, proptest};

                proptest! {
                    #[test]
                    fn ldl(m in dmatrix_($scalar)) {
                        let m = &m * m.adjoint();

                        if let Some(ldl) = m.clone().ldl() {
                            let p = &ldl.l * &ldl.d_matrix() * &ldl.l.transpose();
                            println!("m: {}, p: {}", m, p);

                            prop_assert!(relative_eq!(m, p, epsilon = 1.0e-7));
                        }
                    }

                    #[test]
                    fn ldl_static(m in matrix4_($scalar)) {
                        let m = m.hermitian_part();

                        if let Some(ldl) = m.ldl() {
                            let p = ldl.l * ldl.d_matrix() * ldl.l.transpose();
                            prop_assert!(relative_eq!(m, p, epsilon = 1.0e-7));
                        }
                    }
                }
            }
        }
    );

    gen_tests!(f64, PROPTEST_F64);
}