rusoda 0.1.0

Pure rust and just need log dependencies implementation of odepack dlsoda.
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

use super::{daxpy::daxpy, ddot::ddot};

/// Solves the double precision system A * x = b  or  trans(A) * x = b using the factors computed by dgeco or dgefa.
#[no_mangle]
#[inline(always)]
pub fn dgesl(a: &mut Vec<Vec<f64>>, n: usize, ipvt: &[usize], b: &mut [f64], job: usize) {
    /*
    [2023-10-17]

    Impl of linpack dgesl.f by Elliot Law.

    [Purpose]
    DGESL solves the double precision system
    A * x = b  or  trans(A) * x = b
    using the factors computed by dgeco or dgefa.

    [Author]
    Cleve Moler,
    University of New Mexico,
    Argonne National Lab.

    */
    let mut t: f64;
    let mut j: usize;
    /*
       Job = 0, solve a * x = b.
    */
    if job == 0 {
        /*
           First solve L * y = b.
        */
        for k in 1..n + 1 {
            t = ddot(k - 1, &a[k], b, 1, 1);
            b[k] = (b[k] - t) / a[k][k];
        }
        /*
           Now solve U * x = y.
        */
        for k in (1..n).rev() {
            (*b)[k] += ddot(n - k, &a[k][k..], &b[k..], 1, 1);
            j = ipvt[k];
            if j != k {
                t = b[j];
                b[j] = b[k];
                b[k] = t;
            }
        }
        return;
    }
    /*
       Job = nonzero, solve Transpose(a) * x = b.

       First solve Transpose(U) * y = b.
    */
    for k in 1..n {
        j = ipvt[k];
        t = b[j];
        if j != k {
            b[j] = b[k];
            b[k] = t;
        }
        daxpy(n - k, &mut a[k][k..], &mut b[k..], 1, 1, &t);
    }
    /*
       Now solve Transpose(L) * x = y.
    */
    for k in (1..n + 1).rev() {
        b[k] /= a[k][k];
        t = -b[k];
        daxpy(k - 1, &mut a[k], b, 1, 1, &t);
    }
}