sparsemat 0.2.0

A simple sparse matrix library
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

use crate::types::*;
use crate::sparsematrix::*;
use crate::vector::*;
use crate::densevec::*;

pub trait LinearSolver {
    fn solve<'a, M>(&self, mat: &'a M, b: &DenseVec<M::Value>, x: &mut DenseVec<M::Value>)
    where M: SparseMatrix<'a>,
          M::Value: FloatType;
}

pub struct ConjugateGradient {
    tol: f64,
    iter_max: usize,
}

impl Default for ConjugateGradient {
    fn default() -> Self {
        Self {
            tol: 1e-12f64,
            iter_max: 10_000,
        }
    }
}

impl LinearSolver for ConjugateGradient {
    fn solve<'a, M>(&self, mat: &'a M, b: &DenseVec<M::Value>, x: &mut DenseVec<M::Value>)
    where M: SparseMatrix<'a>,
          M::Value: FloatType {
        if mat.n_rows() != mat.n_cols() {
            panic!("Matrix is not symmetric");
        }
        if mat.n_rows() != b.dim()
            || mat.n_rows() != x.dim() {
            panic!("Matrix and vector size mismatch");
        }
        // r = b - M * x
        let mut r = b.clone() - mat.mvp(x);
        let mut p = r.clone();
        let mut r_norm_squared = r.norm_squared();
        for _k in 0..self.iter_max {
            // M * p
            let mat_p = mat.mvp(&p);
            // alpha = r * r / (p * M * p)
            let alpha = r_norm_squared / p.inner_prod(&mat_p);
            // x = x + alpha * p
            *x += p.clone() * alpha;
            // r = r - alpha * (M * p)
            r -= mat_p * alpha;
            let r_norm_squared_prev = r_norm_squared;
            r_norm_squared = r.norm_squared();
            if f64::sqrt(r_norm_squared.into()) < self.tol {
                break;
            }
            // beta = r * r / (r_prev * r_prev)
            let beta = r_norm_squared / r_norm_squared_prev;
            // p = r + beta + p
            p.scale(beta);
            p.add(&r);
        }
    }
}