scimath 0.1.2

A scientific computing library. WIP
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

use crate::types::matrix::Matrix;

// solving system Ax = b by LU factorization
pub fn forward_sub(L: &Matrix, b: &Matrix) -> Matrix {
    let n: usize = L.shape().0;
    let mut x = Matrix::new(n, 1);
    let mut b_: Matrix = b.clone();
    for j in 0..n {
        if L[(j,j)] == 0.0 {
            panic!("Lower triangular matrix is singular!");
        }
        x[(j,0)] = b_[(j,0)] / L[(0,0)];
        for i in j + 1 .. n {
            b_[(i,0)] = b_[(i,0)] - L[(i,j)] * x[(j,0)];
        }
    }
    return x;
}

pub fn back_sub(U: &Matrix, b: &Matrix) -> Matrix {
    let n: usize = U.shape().0;
    let mut x = Matrix::new(n, 1);
    let mut b_ = b.clone();

    for j in (0..n).rev() {
        if U[(j,j)] == 0.0 {
            panic!("Upper triangular matrix is singular!");
        }
        x[(j,0)] = b_[(j,0)] / U[(j,j)];
        for i in 0..j {
            b_[(i,0)] = b_[(i,0)] - U[(i,j)] * x[(j,0)];
        }
    }
    return x;
}

pub fn lu(A: &Matrix, b:&Matrix) -> Matrix {
    let n = A.shape().0;
    let mut LU = A.clone();

    for k in 0..n {
        if LU[(k,k)] == 0.0 {
            panic!("LU factorization does not exist!");
        }
        for i in k + 1 .. n {
            LU[(i,k)] = LU[(i,k)] / LU[(k,k)];
        }
        for j in k + 1 .. n {
            for i in k + 1 .. n {
                LU[(i,j)] = LU[(i,j)] - LU[(i,k)] * LU[(k,j)];
            }
        }
    }

    return LU;
}

pub fn partial_lu(A: &Matrix, b:&Matrix) -> (Matrix, Matrix) {
    let n = A.shape().0;
    let mut LU = A.clone();
    let mut Pb = b.clone();

    for k in 0..n {
        // find p s.t. pivot
        let mut p: usize = 0;
        let mut max: f64 = 0.0;
        for i in k..n {
            if LU[(i,k)].abs() > max {
                p = i;
                max = LU[(i,k)].abs();
            }
        }
        if p != k {
            LU.swap_row(p,k);
            Pb.swap_row(p,k);
        }
        if LU[(k,k)] == 0.0 {
            continue;
        }
        for i in k + 1 .. n {
            LU[(i,k)] = LU[(i,k)] / LU[(k,k)];
        }
        for j in k + 1 .. n {
            for i in k + 1 .. n {
                LU[(i,j)] = LU[(i,j)] - LU[(i,k)] * LU[(k,j)];
                Pb[(i,0)] = Pb[(i,0)] - LU[(i,k)] * LU[(k,j)];
            }
        }
    }

    return (LU,Pb);
}

pub fn solve(A: &Matrix, b: &Matrix) -> Matrix {
    let LU = partial_lu(A, b);
    let x = back_sub(&LU.0, &LU.1);
    return x;
}