ldl 0.1.0

LDL factorisation for quasi-definite linear systems
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

extern crate ldl;

fn main() {
    let a_n: usize = 10;
    let a_p = vec![0, 1, 2, 4, 5, 6, 8, 10, 12, 14, 17];
    let a_i = vec![0, 1, 1, 2, 3, 4, 1, 5, 0, 6, 3, 7, 6, 8, 1, 2, 9];
    let a_x = vec![
        1.0, 0.460641, -0.121189, 0.417928, 0.177828, 0.1, -0.0290058, -1.0, 0.350321, -0.441092,
        -0.0845395, -0.316228, 0.178663, -0.299077, 0.182452, -1.56506, -0.1,
    ];
    let b = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0];

    let l_n = a_n;

    // pre-factorisation memory allocations //

    // These can happen *before* the etree is calculated
    // since the sizes are not sparsity pattern specific

    // For the elimination tree
    let mut etree = vec![None; a_n];
    let mut l_nz = vec![0; a_n];

    // For the L factors. Li and Lx are sparsity dependent
    // so must be done after the etree is constructed
    let mut l_p = vec![0; a_n + 1];
    let mut d = vec![0.0; a_n];
    let mut d_inv = vec![0.0; a_n];

    // Working memory. Note that both the etree and factor
    // calls requires a working vector of int, with
    // the factor function requiring 3*An elements and the
    // etree only An elements. Just allocate the larger
    // amount here and use it in both places
    let mut iwork = vec![0; 3 * a_n];
    let mut bwork = vec![ldl::Marker::Unused; a_n];
    let mut fwork = vec![0.0; a_n];

    // Elimination tree calculation //

    let sum_l_nz = ldl::etree::<usize>(a_n, &a_p, &a_i, &mut iwork, &mut l_nz, &mut etree).unwrap();

    // LDL factorisation //

    let mut l_i = vec![0; sum_l_nz as usize];
    let mut l_x = vec![0.0; sum_l_nz as usize];

    ldl::factor::<f64, usize>(
        a_n, &a_p, &a_i, &a_x, &mut l_p, &mut l_i, &mut l_x, &mut d, &mut d_inv, &l_nz, &etree,
        &mut bwork, &mut iwork, &mut fwork,
    )
    .unwrap();

    // Solve //

    let mut x = vec![0.0; a_n as usize];

    // when solving A\b, start with x = b
    for i in 0..l_n {
        x[i as usize] = b[i as usize];
    }
    ldl::solve::<f64, usize>(l_n, &l_p, &l_i, &l_x, &d_inv, &mut x);

    // Print factors and solution //

    println!();
    println!("A (CSC format):");
    println!("--------------------------");
    println!("A.p = {:?}", a_p);
    println!("A.i = {:?}", a_i);
    println!("A.x = {:?}", a_x);
    println!();
    println!();

    println!("elimination tree:");
    println!("--------------------------");
    println!("etree = {:?}", etree);
    println!("Lnz = {:?}", l_nz);
    println!();
    println!();

    println!("L (CSC format):");
    println!("--------------------------");
    println!("L.p = {:?}", l_p);
    println!("L.i = {:?}", l_i);
    println!("L.x = {:?}", l_x);
    println!();
    println!();

    println!("D:");
    println!("--------------------------");
    println!("diag(D)      = {:?}", d);
    println!("diag(D^{{-1}}) = {:?}", d_inv);
    println!();
    println!();

    println!("solve results:");
    println!("--------------------------");
    println!("b = {:?}", b);
    println!("A\\b = {:?}", x);
    println!();
    println!();
}