russell_sparse 2.1.0

Solvers for large sparse linear systems (wraps MUMPS and UMFPACK)
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

#[cfg(feature = "local_sparse")]
use russell_lab::*;

#[cfg(feature = "local_sparse")]
use russell_sparse::prelude::*;

#[cfg(feature = "local_sparse")]
use russell_sparse::StrError;

#[cfg(feature = "local_sparse")]
fn main() -> Result<(), StrError> {
    // constants
    let ndim = 5; // number of rows = number of columns
    let nnz = 13; // number of non-zero values, including duplicates

    // allocate solver
    let mut mumps = SolverMUMPS::new()?;

    // allocate the coefficient matrix
    //  2  3  .  .  .
    //  3  .  4  .  6
    //  . -1 -3  2  .
    //  .  .  1  .  .
    //  .  4  2  .  1
    let mut coo = CooMatrix::new(ndim, ndim, nnz, Sym::No)?;
    coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate
    coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate
    coo.put(1, 0, 3.0)?;
    coo.put(0, 1, 3.0)?;
    coo.put(2, 1, -1.0)?;
    coo.put(4, 1, 4.0)?;
    coo.put(1, 2, 4.0)?;
    coo.put(2, 2, -3.0)?;
    coo.put(3, 2, 1.0)?;
    coo.put(4, 2, 2.0)?;
    coo.put(2, 3, 2.0)?;
    coo.put(1, 4, 6.0)?;
    coo.put(4, 4, 1.0)?;

    // parameters
    let mut params = LinSolParams::new();
    params.compute_error_estimates = true;
    params.compute_condition_numbers = true;
    params.verbose = false;

    // call factorize
    mumps.factorize(&coo, Some(params))?;

    // allocate x and rhs
    let mut x = Vector::new(ndim);
    let rhs = Vector::from(&[8.0, 45.0, -3.0, 3.0, 19.0]);

    // calculate the solution
    mumps.solve(&mut x, &rhs, false)?;
    println!("x =\n{}", x);

    // check the results
    let correct = vec![1.0, 2.0, 3.0, 4.0, 5.0];
    vec_approx_eq(&x, &correct, 1e-14);

    // analysis
    let a = coo.as_dense();
    let mut ai = Matrix::new(5, 5);
    let det_a = mat_inverse(&mut ai, &a).unwrap();
    let norm_a = mat_norm(&a, Norm::Inf);
    let norm_ai = mat_norm(&ai, Norm::Inf);
    let cond = norm_a * norm_ai;
    let rcond = 1.0 / cond;
    let verify = VerifyLinSys::from(&coo, &x, &rhs).unwrap();
    let mut stats = StatsLinSol::new();
    mumps.update_stats(&mut stats);
    let s = &stats.mumps_stats;
    println!("\n___ ANALYSIS ________________________");
    println!("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
    println!("a =\n{}", a);
    println!("a⁻¹ =\n{:.5}", ai);
    println!("              det(a) = {:?}", det_a);
    println!("                ‖a‖∞ = {:?}", norm_a);
    println!("              ‖a⁻¹‖∞ = {:?}", norm_ai);
    println!("cond = ‖a‖∞ · ‖a⁻¹‖∞ = {:?}", cond);
    println!("    rcond = 1 / cond = {:?}", rcond);
    println!("           max_abs_a = {:?}", verify.max_abs_a);
    println!("          max_abs_ax = {:?}", verify.max_abs_ax);
    println!("        max_abs_diff = {:?}", verify.max_abs_diff);
    println!("      relative_error = {:?}", verify.relative_error);
    println!("              norm_a = {:?}", &s.inf_norm_a);
    println!("              norm_x = {:?}", &s.inf_norm_x);
    println!("            residual = {:?}", &s.scaled_residual);
    println!("              omega1 = {:?}", &s.backward_error_omega1);
    println!("              omega2 = {:?}", &s.backward_error_omega2);
    println!("               delta = {:?}", &s.normalized_delta_x);
    println!("               cond1 = {:?}", &s.condition_number1);
    println!("               cond2 = {:?}", &s.condition_number2);
    Ok(())
}

#[cfg(not(feature = "local_sparse"))]
fn main() {
    println!("MUMPS is not available");
}