spsolve 0.1.0

Solve sparse systems of linear equations
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

use crate::{FactorSolver, Solver};
use anyhow::{format_err, Result};
use num_traits::{NumAssignOps, PrimInt};
use rlu::{Scalar, LU};
use std::fmt::Display;

/// Solver based on [AMD](https://crates.io/crates/amd) and [RLU](https://crates.io/crates/rlu).
#[derive(Default)]
pub struct RLU {
    pub control: amd::Control,
    pub options: rlu::Options,
}

impl<I, S> Solver<I, S> for RLU
where
    I: PrimInt + NumAssignOps + Display,
    S: Scalar,
{
    fn solve(
        &self,
        n: usize,
        a_i: &[I],
        a_p: &[I],
        a_x: &[S],
        b: &mut [S],
        trans: bool,
    ) -> Result<()> {
        let (p, _p_inv, _info) = amd::order::<I>(I::from(n).unwrap(), &a_p, &a_i, &self.control)
            .map_err(|st| format_err!("amd status: {:?}", st))?;

        let lu = rlu::factor(
            I::from(n).unwrap(),
            &a_i,
            &a_p,
            &a_x,
            Some(&p),
            &self.options,
        )
        .map_err(|err| format_err!("factor error: {}", err))?;

        rlu::solve(&lu, b, trans).map_err(|err| format_err!("solve error: {}", err))?;

        Ok(())
    }
}

impl<I, S> FactorSolver<I, S, LU<S>> for RLU
where
    I: PrimInt + NumAssignOps + Display,
    S: Scalar,
{
    fn factor(&self, n: usize, a_i: &[I], a_p: &[I], a_x: &[S]) -> Result<LU<S>> {
        let (p, _p_inv, _info) = amd::order::<I>(I::from(n).unwrap(), &a_p, &a_i, &self.control)
            .map_err(|st| format_err!("amd status: {:?}", st))?;

        let lu = rlu::factor(
            I::from(n).unwrap(),
            &a_i,
            &a_p,
            &a_x,
            Some(&p),
            &self.options,
        )
        .map_err(|err| format_err!("factor error: {}", err))?;
        Ok(lu)
    }

    fn solve(&self, f: &LU<S>, b: &mut [S], trans: bool) -> Result<()> {
        rlu::solve(f, b, trans).map_err(|err| format_err!("solve error: {}", err))?;
        Ok(())
    }
}

#[cfg(test)]
mod tests {
    use super::RLU;
    use crate::test;
    use anyhow::Result;

    #[test]
    fn simple_test() -> Result<()> {
        let solver = RLU::default();
        test::simple_solver_test::<usize, f64, RLU>(solver)
    }

    #[test]
    fn test_solver() -> Result<()> {
        let solver = RLU::default();
        test::test_bbus(&solver, 1, 1e-11)
    }
}