From f859b26edbfc8102fabc67b8a417bfcfa9c6ace5 Mon Sep 17 00:00:00 2001 From: laentropia Date: Wed, 29 Apr 2026 16:53:28 -0600 Subject: [PATCH] test: gaussian and gauus jordan elimination --- src/lib.rs | 446 +++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 436 insertions(+), 10 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 9896c49..3a88683 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,4 +1,5 @@ use fractions::{Fraction, FractionError}; +use std::cmp::min; use std::fmt::{Debug, Display}; use std::ops::Add; use std::ops::Mul; @@ -232,10 +233,10 @@ impl Matrix { } let mut max_row = col; - let mut max_value = self.get(col, col).unwrap().abs(); + let mut max_value = self.get(col, col)?.abs(); for r in (col + 1)..self.rows { - let val = self.get(r, col).unwrap().abs(); + let val = self.get(r, col)?.abs(); if val > max_value { max_value = val; max_row = r; @@ -260,24 +261,41 @@ impl Matrix { rows: self.rows, data: self.data.clone(), }; - let mut sign = Fraction::new(1, 1).unwrap(); + let mut sign = Fraction::new(1, 1)?; - for i in 0..self.columns { + for i in 0..min(self.rows, self.columns) { // We do partial pivoting for better efifiency and security let pivot_exists = trig_matrix.partial_pivoting(i, &mut sign)?; // If there ain't no other thing but 0 then we're // fucked, determinant is zero if !pivot_exists { - return Err(MatrixError::FailedGauss); + let mut all_zero = true; + + for r in i..trig_matrix.rows { + if !trig_matrix.get(r, i)?.is_zero() { + all_zero = false; + break; + } + } + + if all_zero { + return Err(MatrixError::FailedGauss); // matriz singular + } + + continue; } - let pivot = *trig_matrix.get(i, i).unwrap(); + let pivot = *trig_matrix.get(i, i)?; + + if pivot.is_zero() { + return Err(MatrixError::FailedGauss); + } // The main gaussian elimination, not even I remember how // i did it in such a asimple way for x in (i + 1)..trig_matrix.rows { - let m = (*trig_matrix.get(x, i).unwrap() / pivot).unwrap(); + let m = (*trig_matrix.get(x, i)? / pivot)?; let row_x = trig_matrix.get_row(x)?; let row_i = trig_matrix.get_row(i)?; @@ -303,14 +321,14 @@ impl Matrix { let mut dummy = Fraction::from(1); - for i in 0..self.columns { + for i in 0..min(self.columns, self.rows) { let pivot_exists = new_matrix.partial_pivoting(i, &mut dummy)?; if !pivot_exists { return Err(MatrixError::FailedGaussJordan); } - let pivot = *new_matrix.get(i, i).unwrap(); + let pivot = *new_matrix.get(i, i)?; let new_pivot_row = new_matrix .get_row(i)? @@ -324,7 +342,7 @@ impl Matrix { continue; } - let factor = *new_matrix.get(r, i).unwrap(); + let factor = *new_matrix.get(r, i)?; if factor.is_zero() { continue; @@ -2128,4 +2146,412 @@ mod tests { } } } + + #[test] + fn test_partial_pivoting_no_swap_needed() { + let mut m = Matrix { + rows: 3, + columns: 3, + data: vec![ + Fraction::from(5), + Fraction::from(0), + Fraction::from(0), + Fraction::from(2), + Fraction::from(0), + Fraction::from(0), + Fraction::from(1), + Fraction::from(0), + Fraction::from(0), + ], + }; + + let mut sign = Fraction::from(1); + + let res = m.partial_pivoting(0, &mut sign).unwrap(); + + assert!(res); + assert_eq!(sign, Fraction::from(1)); // no swap + } + + #[test] + fn test_partial_pivoting_with_swap() { + let mut m = Matrix { + rows: 3, + columns: 3, + data: vec![ + Fraction::from(1), + Fraction::from(0), + Fraction::from(0), + Fraction::from(9), + Fraction::from(0), + Fraction::from(0), + Fraction::from(3), + Fraction::from(0), + Fraction::from(0), + ], + }; + + let mut sign = Fraction::from(1); + + let res = m.partial_pivoting(0, &mut sign).unwrap(); + + assert!(res); + + // ahora la fila 0 debe ser la de valor 9 + assert_eq!(m.get(0, 0).unwrap(), &Fraction::from(9)); + + assert_eq!(sign, Fraction::from(-1)); // hubo swap + } + + #[test] + fn test_partial_pivoting_all_zero_column() { + let mut m = Matrix { + rows: 3, + columns: 3, + data: vec![ + Fraction::from(0), + Fraction::from(1), + Fraction::from(2), + Fraction::from(0), + Fraction::from(3), + Fraction::from(4), + Fraction::from(0), + Fraction::from(5), + Fraction::from(6), + ], + }; + + let mut sign = Fraction::from(1); + + let res = m.partial_pivoting(0, &mut sign).unwrap(); + + assert!(!res); // no pivot posible + } + + #[test] + fn test_partial_pivoting_largest_abs_value() { + let mut m = Matrix { + rows: 3, + columns: 3, + data: vec![ + Fraction::from(-3), + Fraction::from(0), + Fraction::from(0), + Fraction::from(7), + Fraction::from(0), + Fraction::from(0), + Fraction::from(-10), + Fraction::from(0), + Fraction::from(0), + ], + }; + + let mut sign = Fraction::from(1); + + m.partial_pivoting(0, &mut sign).unwrap(); + + assert_eq!( + m.get(0, 0).unwrap(), + &Fraction::from(-10) // mayor en valor absoluto + ); + } + + #[test] + fn test_gaussian_already_triangular() { + let m = Matrix { + rows: 3, + columns: 3, + data: vec![ + Fraction::from(1), + Fraction::from(2), + Fraction::from(3), + Fraction::from(0), + Fraction::from(4), + Fraction::from(5), + Fraction::from(0), + Fraction::from(0), + Fraction::from(6), + ], + }; + + let (res, _) = m.gaussian_elimination().unwrap(); + + assert_eq!(res.data, m.data); + } + + #[test] + fn test_gaussian_simple() { + let m = Matrix { + rows: 2, + columns: 2, + data: vec![ + Fraction::from(2), + Fraction::from(1), + Fraction::from(4), + Fraction::from(3), + ], + }; + + let (res, _) = m.gaussian_elimination().unwrap(); + + // abajo izquierda debe ser 0 + assert!(res.get(1, 0).unwrap().is_zero()); + } + + #[test] + fn test_gaussian_requires_pivoting() { + let m = Matrix { + rows: 2, + columns: 2, + data: vec![ + Fraction::from(0), + Fraction::from(1), + Fraction::from(2), + Fraction::from(3), + ], + }; + + let (res, _) = m.gaussian_elimination().unwrap(); + + assert!(res.get(0, 0).unwrap() != &Fraction::from(0)); + } + + #[test] + fn test_gaussian_singular_matrix() { + let m = Matrix { + rows: 2, + columns: 2, + data: vec![ + Fraction::from(1), + Fraction::from(2), + Fraction::from(2), + Fraction::from(4), + ], + }; + + let res = m.gaussian_elimination(); + + assert!(matches!(res, Err(MatrixError::FailedGauss))); + } + + #[test] + fn test_gaussian_upper_triangular_property() { + let m = Matrix { + rows: 3, + columns: 3, + data: vec![ + Fraction::from(2), + Fraction::from(1), + Fraction::from(3), + Fraction::from(4), + Fraction::from(1), + Fraction::from(6), + Fraction::from(6), + Fraction::from(2), + Fraction::from(9), + ], + }; + + let res = m.gaussian_elimination(); + + assert!(res.is_err()); + } + + #[test] + fn test_gj_identity_matrix() { + let m = Matrix { + rows: 3, + columns: 3, + data: vec![ + Fraction::from(1), + Fraction::from(0), + Fraction::from(0), + Fraction::from(0), + Fraction::from(1), + Fraction::from(0), + Fraction::from(0), + Fraction::from(0), + Fraction::from(1), + ], + }; + + let res = m.gauss_jordan_elimination().unwrap(); + + assert_eq!(res.data, m.data); + } + + #[test] + fn test_gj_already_reduced() { + let m = Matrix { + rows: 2, + columns: 2, + data: vec![ + Fraction::from(1), + Fraction::from(0), + Fraction::from(0), + Fraction::from(1), + ], + }; + + let res = m.gauss_jordan_elimination().unwrap(); + + assert_eq!(res.data, m.data); + } + + #[test] + fn test_gj_simple() { + let m = Matrix { + rows: 2, + columns: 2, + data: vec![ + Fraction::from(2), + Fraction::from(1), + Fraction::from(4), + Fraction::from(3), + ], + }; + + let res = m.gauss_jordan_elimination().unwrap(); + + // debe quedar identidad + assert_eq!( + res.data, + vec![ + Fraction::from(1), + Fraction::from(0), + Fraction::from(0), + Fraction::from(1), + ] + ); + } + + #[test] + fn test_gj_requires_pivoting() { + let m = Matrix { + rows: 2, + columns: 2, + data: vec![ + Fraction::from(0), + Fraction::from(1), + Fraction::from(2), + Fraction::from(3), + ], + }; + + let res = m.gauss_jordan_elimination().unwrap(); + + // diagonal debe ser 1 + assert_eq!(*res.get(0, 0).unwrap(), Fraction::from(1)); + assert_eq!(*res.get(1, 1).unwrap(), Fraction::from(1)); + + // fuera de diagonal debe ser 0 + assert!(res.get(0, 1).unwrap().is_zero()); + assert!(res.get(1, 0).unwrap().is_zero()); + } + + #[test] + fn test_gj_3x3() { + let m = Matrix { + rows: 3, + columns: 3, + data: vec![ + Fraction::from(1), + Fraction::from(2), + Fraction::from(3), + Fraction::from(0), + Fraction::from(1), + Fraction::from(4), + Fraction::from(5), + Fraction::from(6), + Fraction::from(0), + ], + }; + + let res = m.gauss_jordan_elimination().unwrap(); + + // debe quedar identidad + for i in 0..3 { + for j in 0..3 { + if i == j { + assert_eq!(*res.get(i, j).unwrap(), Fraction::from(1)); + } else { + assert!(res.get(i, j).unwrap().is_zero()); + } + } + } + } + + #[test] + fn test_gj_singular_matrix() { + let m = Matrix { + rows: 2, + columns: 2, + data: vec![ + Fraction::from(1), + Fraction::from(2), + Fraction::from(2), + Fraction::from(4), + ], + }; + + let res = m.gauss_jordan_elimination(); + + assert!(matches!(res, Err(MatrixError::FailedGaussJordan))); + } + + #[test] + fn test_gj_reduced_row_echelon_form() { + let m = Matrix { + rows: 3, + columns: 3, + data: vec![ + Fraction::from(2), + Fraction::from(1), + Fraction::from(-1), + Fraction::from(-3), + Fraction::from(-1), + Fraction::from(2), + Fraction::from(-2), + Fraction::from(1), + Fraction::from(2), + ], + }; + + let res = m.gauss_jordan_elimination().unwrap(); + + for i in 0..3 { + // pivote = 1 + assert_eq!(*res.get(i, i).unwrap(), Fraction::from(1)); + + // columna del pivote = 0 excepto en pivote + for r in 0..3 { + if r != i { + assert!(res.get(r, i).unwrap().is_zero()); + } + } + } + } + + #[test] + fn test_gj_rectangular_matrix() { + let m = Matrix { + rows: 2, + columns: 3, + data: vec![ + Fraction::from(1), + Fraction::from(2), + Fraction::from(3), + Fraction::from(4), + Fraction::from(5), + Fraction::from(6), + ], + }; + + let res = m.gauss_jordan_elimination().unwrap(); + + // al menos debe cumplir forma escalonada reducida parcial + for i in 0..2 { + assert_eq!(*res.get(i, i).unwrap(), Fraction::from(1)); + } + } }