la-stack 0.4.6

Fast, stack-allocated linear algebra for fixed dimensions
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
142
143
144
145
146

#![forbid(unsafe_code)]

//! Regression coverage for finite solve results and substitution error order.

use la_stack::{ArithmeticOperation, DEFAULT_SINGULAR_TOL, LaError, Matrix, Vector};
use pastey::paste;

fn assert_diagonal_overflow<const D: usize>() {
    // Exercise every output coordinate, including the last one finalized by LU.
    for index in 0..D {
        let mut rows = Matrix::<D>::identity().into_rows();
        rows[index][index] = 1.0e-11;
        let a = Matrix::try_from_rows(rows).unwrap();
        let mut rhs = [0.0; D];
        rhs[index] = 1.0e300;
        let b = Vector::try_new(rhs).unwrap();
        assert_eq!(
            a.lu(DEFAULT_SINGULAR_TOL).unwrap().solve(b),
            Err(LaError::non_finite_computation_step(
                ArithmeticOperation::LuSolve,
                index,
            )),
        );
        assert_eq!(
            a.ldlt(DEFAULT_SINGULAR_TOL).unwrap().solve(b),
            Err(LaError::non_finite_computation_step(
                ArithmeticOperation::LdltSolve,
                index,
            )),
        );
    }
}

fn assert_forward_overflow<const D: usize>() {
    let mut rows = Matrix::<D>::identity().into_rows();
    rows[D - 1][0] = -1.0;
    let mut rhs = [0.0; D];
    rhs[0] = f64::MAX;
    rhs[D - 1] = f64::MAX;
    let b = Vector::try_new(rhs).unwrap();
    let lu = Matrix::try_from_rows(rows)
        .unwrap()
        .lu(DEFAULT_SINGULAR_TOL)
        .unwrap();
    assert_eq!(
        lu.solve(b),
        Err(LaError::non_finite_computation_step(
            ArithmeticOperation::LuSolve,
            D - 1,
        )),
    );

    // A = L Lᵀ with L[D-1, 0] = -1 and unit diagonal, so the same
    // forward substitution overflows before either solve reaches finalization.
    rows[0][D - 1] = -1.0;
    rows[D - 1][D - 1] = 2.0;
    let ldlt = Matrix::try_from_rows(rows)
        .unwrap()
        .ldlt(DEFAULT_SINGULAR_TOL)
        .unwrap();
    assert_eq!(
        ldlt.solve(b),
        Err(LaError::non_finite_computation_step(
            ArithmeticOperation::LdltSolve,
            D - 1,
        )),
    );
}

fn assert_back_substitution_overflow<const D: usize>() {
    // Test each possible failing row. In the all-rows case, a final ascending
    // scan would report 0 instead of the required descending solve step D-2.
    for failing_row in 0..D - 1 {
        for all_rows in [false, true] {
            let mut upper = Matrix::<D>::identity().into_rows();
            let mut spd = upper;
            let mut last_diagonal = 1.0;
            for i in 0..D - 1 {
                if all_rows || i == failing_row {
                    upper[i][D - 1] = 2.0;
                    spd[i][D - 1] = 2.0;
                    spd[D - 1][i] = 2.0;
                    last_diagonal += 4.0;
                }
            }
            // A = L Lᵀ, where L's final row contains the selected multipliers.
            spd[D - 1][D - 1] = last_diagonal;
            let mut rhs = [0.0; D];
            rhs[D - 1] = 1.0e308;
            let b = Vector::try_new(rhs).unwrap();
            let expected_index = if all_rows { D - 2 } else { failing_row };
            let lu = Matrix::try_from_rows(upper)
                .unwrap()
                .lu(DEFAULT_SINGULAR_TOL)
                .unwrap();
            assert_eq!(
                lu.solve(b),
                Err(LaError::non_finite_computation_step(
                    ArithmeticOperation::LuSolve,
                    expected_index,
                )),
            );
            let ldlt = Matrix::try_from_rows(spd)
                .unwrap()
                .ldlt(DEFAULT_SINGULAR_TOL)
                .unwrap();
            assert_eq!(
                ldlt.solve(b),
                Err(LaError::non_finite_computation_step(
                    ArithmeticOperation::LdltSolve,
                    expected_index,
                )),
            );
        }
    }
}

macro_rules! gen_tests {
    ($d:literal) => {
        paste! {
            #[test]
            fn [<diagonal_overflow_ $d d>]() {
                assert_diagonal_overflow::<$d>();
            }

            #[test]
            fn [<forward_overflow_ $d d>]() {
                assert_forward_overflow::<$d>();
            }

            #[test]
            fn [<back_substitution_overflow_ $d d>]() {
                assert_back_substitution_overflow::<$d>();
            }
        }
    };
}

gen_tests!(2);
gen_tests!(3);
gen_tests!(4);
gen_tests!(5);
gen_tests!(8);
gen_tests!(16);
gen_tests!(32);
gen_tests!(64);