rexl_matrix/vec/
vec_matrix.rs

1/// row based
2pub struct VecMatrix<T> {
3    pub rows: usize,
4    pub cols: usize,
5    pub data: Vec<Vec<T>>,
6}
7
8impl VecMatrix<f64> {
9    pub fn new(rows: usize, cols: usize) -> VecMatrix<f64> {
10        VecMatrix { rows, cols, data: vec![vec![0.; cols]; rows] }
11    }
12
13    pub fn pascal(level: usize) -> VecMatrix<f64> {
14        let mut this = VecMatrix::new(level, level);
15        for i in 0..level {
16            for j in 0..level {
17                if i == 0 || j == 0 {
18                    this.data[i][j] = 1.
19                } else {
20                    this.data[i][j] = this.data[i - 1][j] + this.data[i][j - 1];
21                }
22            }
23        }
24        this
25    }
26
27    pub fn swap_row(&mut self, i1: usize, i2: usize) {
28        self.data.swap(i1, i2);
29    }
30
31    pub fn add_other_row(&mut self, dest: usize, src: usize, multiply: f64) {
32        for i in 0..self.cols {
33            self.data[dest][i] += self.data[src][i] * multiply;
34        }
35    }
36
37    pub fn multiply_row(&mut self, dest: usize, multiply: f64) {
38        for i in 0..self.cols {
39            self.data[dest][i] *= multiply;
40        }
41    }
42
43    pub fn diagonalize(&mut self) {
44        let mut dimension = self.rows;
45        if self.rows > self.cols {
46            dimension = self.cols
47        }
48
49        let mut sign = true;
50        for k in 0..dimension - 1 {
51            if self.data[k][k] == 0. {
52                let mut swapped = false;
53                // swap rows to make it no-zero
54                for i in (k + 1)..dimension {
55                    if self.data[i][k] != 0. {
56                        self.swap_row(k, i);
57                        swapped = true;
58                        // reverse sign
59                        sign = !sign;
60                        break;
61                    }
62                }
63                // cannot be diagonalized
64                if !swapped {
65                    return;
66                }
67            }
68
69            for i in (k + 1)..dimension {
70                self.add_other_row(i, k, -self.data[i][k] / self.data[k][k]);
71            }
72        }
73
74        if !sign {
75            self.multiply_row(0, -1.);
76        }
77    }
78
79    pub fn det(&mut self) -> f64 {
80        let mut dimension = self.rows;
81        if self.rows > self.cols {
82            dimension = self.cols
83        }
84
85        self.diagonalize();
86        let mut result = 1.;
87        for i in 0..dimension {
88            result *= self.data[i][i];
89        }
90        result
91    }
92
93    pub fn print_data(&self) {
94        for row in &self.data {
95            println!("{:?}", row);
96        }
97    }
98}