math_rs/matrix/float/
inverse.rs

1use crate::matrix::traits::{Identity, Invertible, Matrix};
2
3use super::MatrixF32;
4
5impl Invertible for MatrixF32 {
6    fn inverse_gauss_jordan(&self) -> crate::Result<Self>
7    where
8        Self: Sized,
9    {
10        let mut matrix = self.clone(); // I don't want to lose the original matrix
11        let mut inverse = MatrixF32::id(self.rows(), self.tolerance());
12        for i in 0..self.rows() {
13            // Find the pivot
14            let mut pivot_row = i;
15            for j in i + 1..self.rows() {
16                if matrix.get(j, i)?.abs() > matrix.get(pivot_row, i)?.abs() {
17                    pivot_row = j;
18                }
19            }
20
21            // Swap the rows
22            if pivot_row != i {
23                matrix.swap_rows(i, pivot_row)?;
24                inverse.swap_rows(i, pivot_row)?;
25            }
26
27            // Divide the row by the pivot
28            let pivot = matrix.get(i, i)?.clone();
29            if pivot.abs() < self.tolerance() {
30                return Err(crate::MathError::MatrixError(
31                    "Matrix is not invertible".to_string(),
32                ));
33            }
34            for j in 0..self.columns() {
35                matrix.set(i, j, matrix.get(i, j)? / pivot)?;
36                inverse.set(i, j, inverse.get(i, j)? / pivot)?;
37            }
38
39            // Subtract the row from all other rows
40            for j in 0..self.rows() {
41                if j != i {
42                    let factor = matrix.get(j, i)?.to_owned();
43                    for k in 0..self.columns() {
44                        let new_value = matrix.get(j, k)? - factor * matrix.get(i, k)?;
45                        matrix.set(j, k, new_value)?;
46                        let new_value = inverse.get(j, k)? - factor * inverse.get(i, k)?;
47                        inverse.set(j, k, new_value)?;
48                    }
49                }
50            }
51        }
52        Ok(inverse)
53    }
54
55    fn inverse_montante(&self) -> crate::Result<Self>
56    where
57        Self: Sized,
58    {
59        todo!()
60    }
61
62    fn inverse_adjoint(&self) -> crate::Result<Self>
63    where
64        Self: Sized,
65    {
66        todo!()
67    }
68}
69
70#[cfg(test)]
71mod test {
72    use crate::matrix::traits::{Invertible, Parseable};
73    use crate::matrix::MatrixF32;
74    use crate::matrix_f32;
75
76    const TOL: f32 = 1e-4;
77
78    #[test]
79    fn inverse_gauss_jordan_2x2_f32() {
80        let mat_a = matrix_f32!("{{1,2},{3,4}}", TOL).unwrap();
81        let mat_a_inv = mat_a.inverse_gauss_jordan().unwrap();
82        let mat_a_inv_expected = matrix_f32!("{{-2,1},{1.5,-0.5}}", TOL).unwrap();
83        pretty_assertions::assert_eq!(mat_a_inv, mat_a_inv_expected);
84    }
85
86    #[test]
87    fn inverse_gauss_jordan_3x3_f32() {
88        let mat_a = matrix_f32!("{{1,2,3},{0,1,4},{5,6,0}}", TOL).unwrap();
89        let mat_a_inv = mat_a.inverse_gauss_jordan().unwrap();
90        let mat_a_inv_expected = matrix_f32!("{{-24,18,5},{20,-15,-4},{-5,4,1}}", TOL).unwrap();
91        pretty_assertions::assert_eq!(mat_a_inv, mat_a_inv_expected);
92    }
93
94    #[test]
95    fn inverse_gauss_jordan_4x4_f32() {
96        let mat_a = matrix_f32!("{{1,2,3,4},{0,1,4,5},{5,6,0,7},{8,9,10,0}}", TOL).unwrap();
97        let mat_a_inv = mat_a.inverse_gauss_jordan().unwrap();
98        let mat_a_inv_expected = matrix_f32!(
99            "{{-5.0736837,3.0421052,0.72631574,0.30526316},{4.9894733,-3.1368423,-0.6105262,-0.24210525},{ -0.43157893,0.38947368,-0.03157895,0.07368421},{-0.6526315,0.51578945,0.1473684,-0.010526315}}",
100            TOL
101        )
102        .unwrap();
103        pretty_assertions::assert_eq!(mat_a_inv, mat_a_inv_expected);
104    }
105
106    #[test]
107    fn inverse_gauss_jordan_5x5_f32() {
108        let mat_a = matrix_f32!(
109            "{{1,2,3,4,5},{0,1,4,5,6},{5,6,0,7,8},{8,9,10,0,11},{12,13,14,15,0}}",
110            TOL
111        )
112        .unwrap();
113        let mat_a_inv = mat_a.inverse_gauss_jordan().unwrap();
114        let mat_a_inv_expected = matrix_f32!(
115            "{{-4.25986, +2.42993, +0.55141, +0.20986, +0.06866} ,{+4.42606, -2.71303, -0.48944, -0.17606, -0.04754} ,{-0.13944, +0.16972, -0.09437, +0.03944, +0.02465} ,{-0.29789, +0.24894, +0.07113, -0.05211, +0.02993} ,{-0.39648, +0.29824, +0.08521, +0.04648, -0.03345}}",
116            TOL
117        )
118        .unwrap();
119        pretty_assertions::assert_eq!(mat_a_inv, mat_a_inv_expected);
120    }
121}
math_rs/matrix/float/inverse.rs

math_rs/matrix/float/
inverse.rs
