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

#![forbid(unsafe_code)]

//! Exact linear system solve for a near-singular 3×3 system.
//!
//! This example demonstrates `solve_exact()` and [`ExactF64Conversion`]. The
//! solve uses arbitrary-precision rational arithmetic to compute a solution
//! exact for the stored binary64 inputs, even when the f64 LU solve is wildly
//! inaccurate. `try_to_f64()` only succeeds when every exact component is
//! exactly representable as `f64`, while
//! `to_rounded_f64()` explicitly opts into nearest-even rounding.
//!
//! Run with: `cargo run --features exact --example exact_solve_3x3`

use la_stack::prelude::*;

fn main() -> Result<(), LaError> {
    // Near-singular 3×3 system.
    //
    // The base matrix [[1,2,3],[4,5,6],[7,8,9]] is exactly singular (rows in
    // arithmetic progression).  Perturbing entry (0,0) by 2^-50 ≈ 8.9e-16
    // makes it invertible but extremely ill-conditioned.
    let perturbation = f64::from_bits(0x3CD0_0000_0000_0000); // 2^-50
    let a = Matrix::<3>::try_from_rows([
        [1.0 + perturbation, 2.0, 3.0],
        [4.0, 5.0, 6.0],
        [7.0, 8.0, 9.0],
    ])?;

    let b = Vector::<3>::try_new([1.0, 2.0, 3.0])?;

    // f64 LU solve (using zero pivot tolerance since the matrix is nearly singular
    // and would be rejected by DEFAULT_SINGULAR_TOL).
    let lu_x = a.lu(Tolerance::try_new(0.0)?)?.solve(b)?.into_array();

    // Exact solve.
    let exact_x = a.solve_exact(b)?;
    // Only the third component is nonzero: one third of the third column is b.
    assert_eq!(
        exact_x.as_array(),
        &[
            BigRational::from_integer(0.into()),
            BigRational::from_integer(0.into()),
            BigRational::new(1.into(), 3.into()),
        ],
    );
    println!("Near-singular 3×3 system (perturbation = 2^-50 ≈ {perturbation:.2e}):");
    for row in a.as_rows() {
        print!("  [");
        for (col, value) in row.iter().enumerate() {
            if col > 0 {
                print!(", ");
            }
            print!("{value:22.18}");
        }
        println!("]");
    }
    println!(
        "b = [{}, {}, {}]",
        b.as_array()[0],
        b.as_array()[1],
        b.as_array()[2]
    );
    println!();
    println!(
        "f64 LU solve:       x = [{:+.6e}, {:+.6e}, {:+.6e}]",
        lu_x[0], lu_x[1], lu_x[2]
    );
    println!(
        "solve_exact():      x = [{}, {}, {}]",
        exact_x.as_array()[0],
        exact_x.as_array()[1],
        exact_x.as_array()[2]
    );
    let strict = exact_x.try_to_f64();
    assert_eq!(
        strict,
        Err(LaError::unrepresentable(
            Some(2),
            UnrepresentableReason::RequiresRounding,
        )),
    );
    match strict {
        Ok(x) => {
            let x = x.into_array();
            println!(
                "exact try_to_f64(): x = [{:+.6e}, {:+.6e}, {:+.6e}]",
                x[0], x[1], x[2]
            );
        }
        Err(err) if err.requires_rounding() => {
            println!("exact try_to_f64(): {err}");
            let x = exact_x.to_rounded_f64()?.into_array();
            assert_eq!(x.map(f64::to_bits), [0, 0, 0x3fd5_5555_5555_5555]);
            println!(
                "exact to_rounded_f64(): x = [{:+.6e}, {:+.6e}, {:+.6e}]",
                x[0], x[1], x[2]
            );
        }
        Err(err) => return Err(err),
    }
    Ok(())
}