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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

use nalgebra::{Const, DVector};
use nalgebra_sparse::{io::load_coo_from_matrix_market_file, CsrMatrix};

use crate::{
    io, solver,
    utility::{compute_rel_err, init_b},
};

/// Method is used to set the method used by the routine. 
/// 
/// With:
/// * **Method::JA** = Jacobi
/// * **Method::GS** = Gauss-Seidel
/// * **Method::GR** = gradient
/// * **Method::CG** = conjugate gradient
/// * **Method::PG** = gradient with Jacobi preconditioning
#[derive(Debug)]
pub enum Method {
    JA,
    GS,
    GR,
    CG,
    PG,
}

impl Method {
    pub fn copy(&self) -> Method {
        match self {
            Method::JA => Method::JA,
            Method::GS => Method::GS,
            Method::CG => Method::CG,
            Method::GR => Method::GR,
            Method::PG => Method::PG,
        }
    }
}
/// Performace is a struct returned by [compute_precision].
/// 
/// It contains three fields:
/// * **rel_err**: the relative error of the computed result against the correct one
/// * **time**: the time the routine take to solve the system
/// * **iter**: the number of iteration the routine take to solve the system
pub struct Performance {
    pub rel_err: f64,
    pub time: u128,
    pub iter: u32,
}

/// Stat is a struct returned by [solve_linear_system].
/// 
/// It contains three fields:
/// * **solution**: the solution computed by the routine
/// * **time**: the time the routine take to solve the system
/// * **iter**: the number of iteration the routine take to solve the system
#[derive(Debug)]
pub struct Stat {
    solution: DVector<f64>,
    time: u128,
    iter: u32,
}

impl Stat {
    pub fn new(solution: DVector<f64>, time: u128, iter: u32) -> Stat {
        Stat {
            solution,
            time,
            iter,
        }
    }
    pub fn get_solution(&self) -> &DVector<f64> {
        &self.solution
    }

    pub fn get_time(&self) -> u128 {
        self.time
    }

    pub fn get_iterations(&self) -> u32 {
        self.iter
    }

    pub fn to_string(&self) -> String {
        format!(
            "Result:\n{:?}\nMethod converged in \t{} iterations \t({} ms)",
            self.get_solution(),
            self.get_iterations(),
            self.get_time()
        )
    }
}

/// Parse a matrix market file (.mtx) and return the csr matrix representation of it.
/// Panic if matrix is not sparse
pub fn read_matrix_from_matrix_market_file(file_path: &String) -> CsrMatrix<f64> {
    let coo_matrix = load_coo_from_matrix_market_file(file_path).unwrap();
    CsrMatrix::from(&coo_matrix)
}

/// Parse a vector from a matrix market file or a file with each row representing the entry of the vector.
pub fn read_vector_from_file(file_path: &String) -> DVector<f64> {
    io::parse_vector(file_path)
}

/// Initialize a vector with dimension = \[size\] and each entry = value
pub fn init_solution(size: usize, value: f64) -> DVector<f64> {
    DVector::from_element(size, value)
}

/// Initialize a vector with dimension = \[size\] and each entry a random value between 0 and 1
pub fn init_random_vector(size: usize) -> DVector<f64> {
    DVector::new_random_generic(nalgebra::Dyn(size), Const::<1>)
}

/// Solves the linear system ax=b and returns a [Stat] instance:
/// 
/// ### Arguments:
/// * **a**: matrix
/// * **b**: vector of constant terms
/// * **method**: an istance of the enum Method
/// * **tol**: the tolerance required to stop the routine
/// * **max_iter**: the maximum number of iteration after wich the routine stops
/// * **omega**: the relaxation factor, used only if method is either Jacobi or Gauss-Seidel
pub fn solve_linear_system(
    a: &CsrMatrix<f64>,
    b: &DVector<f64>,
    method: Method,
    tol: f64,
    max_iter: i32,
    omega: f64,
) -> Stat {
    solver::exec(a, b, method, tol, max_iter, omega)
}

/// Determines the accuracy of x (solution computed by the routine) against the correct one and returns an istance of [Performance].
/// 
/// Where x is the solution of the system ax = b, with b := a*solution.
/// 
/// ### Arguments:
/// * **a**: matrix 
/// * **solution**: the given solution of the system
/// * **method**: an istance of the enum Method
/// * **tol**: the tolerance required to stop the routine
/// * **max_iter**: the maximum number of iteration after wich the routine stops
/// * **omega**: the relaxation factor, used only if method is either Jacobi or Gauss-Seidel
pub fn compute_precision(
    a: &CsrMatrix<f64>,
    solution: &DVector<f64>,
    method: Method,
    tol: f64,
    max_iter: i32,
    omega: f64,
) -> Performance {
    let b = init_b(solution, a);
    let result = solver::exec(a, &b, method, tol, max_iter, omega);
    let rel_err = compute_rel_err(solution, result.get_solution());
    Performance {
        rel_err,
        time: result.get_time(),
        iter: result.get_iterations(),
    }
}