solvr 0.2.0

Advanced computing library for real-world problem solving - optimization, differential equations, interpolation, statistics, and more
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

//! Shared sparse iterative solver dispatch.
use crate::DType;

use numr::algorithm::iterative::{
    BiCgStabOptions, CgOptions, GmresOptions, IterativeSolvers, PreconditionerType,
};
use numr::runtime::Runtime;
use numr::sparse::CsrData;
use numr::tensor::Tensor;

use crate::pde::error::{PdeError, PdeResult};
use crate::pde::types::{FdmOptions, Preconditioner, SparseSolver};

/// Result from sparse solve dispatch.
pub struct SparsesolveResult<R: Runtime<DType = DType>> {
    pub solution: Tensor<R>,
    pub iterations: usize,
    pub residual_norm: f64,
}

/// Dispatch to CG/GMRES/BiCGSTAB based on `FdmOptions`, check convergence.
pub fn solve_sparse_system<R, C>(
    client: &C,
    a: &CsrData<R>,
    rhs: &Tensor<R>,
    options: &FdmOptions,
    context: &str,
) -> PdeResult<SparsesolveResult<R>>
where
    R: Runtime<DType = DType>,
    C: IterativeSolvers<R>,
{
    let preconditioner = match options.preconditioner {
        Preconditioner::None => PreconditionerType::None,
        Preconditioner::Ilu0 => PreconditionerType::Ilu0,
        Preconditioner::Ic0 => PreconditionerType::Ic0,
    };

    let (solution, iterations, residual_norm, converged) = match options.solver {
        SparseSolver::Cg => {
            let opts = CgOptions {
                max_iter: options.max_iter,
                rtol: options.tolerance,
                atol: options.tolerance * 0.01,
                preconditioner,
            };
            let r = client.cg(a, rhs, None, opts).map_err(PdeError::from)?;
            (r.solution, r.iterations, r.residual_norm, r.converged)
        }
        SparseSolver::Gmres => {
            let opts = GmresOptions {
                max_iter: options.max_iter,
                rtol: options.tolerance,
                atol: options.tolerance * 0.01,
                restart: 30,
                preconditioner,
                ..Default::default()
            };
            let r = client.gmres(a, rhs, None, opts).map_err(PdeError::from)?;
            (r.solution, r.iterations, r.residual_norm, r.converged)
        }
        SparseSolver::BiCgStab => {
            let opts = BiCgStabOptions {
                max_iter: options.max_iter,
                rtol: options.tolerance,
                atol: options.tolerance * 0.01,
                preconditioner,
            };
            let r = client
                .bicgstab(a, rhs, None, opts)
                .map_err(PdeError::from)?;
            (r.solution, r.iterations, r.residual_norm, r.converged)
        }
    };

    if !converged {
        return Err(PdeError::DidNotConverge {
            iterations,
            tolerance: options.tolerance,
            context: context.to_string(),
        });
    }

    Ok(SparsesolveResult {
        solution,
        iterations,
        residual_norm,
    })
}