kryst 3.2.1

Krylov subspace and preconditioned iterative solvers for dense and sparse linear systems, with shared and distributed memory parallelism.
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

#[cfg(feature = "complex")]
fn main() {
    eprintln!("dense_direct.rs is unavailable when built with --features complex");
}

// Example: Solve a random SPD system using LU and QR direct solvers from kryst.
// to run:
// cargo run --features=dense-direct,mpi --example dense_direct

#[cfg(feature = "dense-direct")]
use faer::Mat;
#[cfg(feature = "dense-direct")]
use kryst::preconditioner::PcSide;
#[cfg(feature = "dense-direct")]
use kryst::solver::legacy::LinearSolver;
#[cfg(feature = "dense-direct")]
use kryst::solver::{LuSolver, QrSolver};
#[cfg(feature = "dense-direct")]
use rand::Rng;

#[cfg(feature = "dense-direct")]
#[cfg(not(feature = "complex"))]
fn main() {
    let comm = kryst::parallel::UniverseComm::NoComm(kryst::parallel::NoComm);
    let n = 10;
    // Build a random SPD matrix: A = MᵀM + I
    let mut rng = rand::thread_rng();
    let data: Vec<f64> = (0..n * n).map(|_| rng.r#gen()).collect();
    let m = Mat::from_fn(n, n, |i, j| data[j * n + i]);
    let m_t = m.transpose();
    // a = m^T * m
    let mut a = &m_t * &m;
    // a = a + I
    for i in 0..n {
        a[(i, i)] = a[(i, i)] + 1.0;
    }

    // Right-hand side
    let b: Vec<f64> = (0..n).map(|_| rng.r#gen()).collect();
    let mut x = vec![0.0; n];

    // LU solve
    let mut lus = LuSolver::new();
    let stats_lu = lus
        .solve(&a, None, &b, &mut x, PcSide::Left, &comm, None, None)
        .unwrap();
    println!("LU x = {:?}, stats = {:?}", x, stats_lu);

    // QR solve
    let mut qrs = QrSolver::new();
    let stats_qr = qrs
        .solve(&a, None, &b, &mut x, PcSide::Left, &comm, None, None)
        .unwrap();
    println!("QR x = {:?}, stats = {:?}", x, stats_qr);
}

#[cfg(not(feature = "dense-direct"))]
#[cfg(not(feature = "complex"))]
fn main() {}