gaussian_elimination.rs

fn gaussian_elimination(a: Vec<Vec<u64>>, b: Vec<u64>) -> Vec<f64> {
    let n = a.len();
    let m = a[0].len();

    let mut augumented: Vec<Vec<f64>> = vec![vec![0.0; m + 1]; n];

    for (y, row) in a.iter().enumerate() {
        for (x, value) in row.iter().enumerate() {
            augumented[y][x] = *value as f64;
        }
    }

    for (y, v) in b.iter().enumerate() {
        augumented[y][m] = *v as f64;
    }

    for i in 0..n {
        let mut pivot_row = i;
        let mut max_val = augumented[i][i].abs();

        for k in i + 1..n {
            if augumented[k][i].abs() > max_val {
                max_val = augumented[k][i];
                pivot_row = k;
            }
        }

        augumented.swap(i, pivot_row);

        let pivot = augumented[i][i];
        if pivot.abs() < 1e-10 {
            panic!("System has no unique solution");
        }

        for v in augumented[i].iter_mut() {
            *v = *v / pivot;
        }

        for j in i + 1..n {
            let factor = augumented[j][i];

            // Split so that 'left' contains everything before index 'j'
            // and 'right' starts exactly at index 'j'.
            let (left, right) = augumented.split_at_mut(j);

            let row_i = &left[i]; // Row i is always in 'left' because i < j
            let row_j = &mut right[0]; // Row j is now the first element of 'right'

            for x in 0..m + 1 {
                row_j[x] -= factor * row_i[x];
            }
        }
    }

    let mut result = vec![0.0; n];

    for i in (0..n).rev() {
        let mut sum = augumented[i][m];

        // Subtract the known values (A_ij * x_j)
        for j in i + 1..n {
            sum -= augumented[i][j] * result[j];
        }

        result[i] = sum;
    }

    return result;
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_gauss_1() {
        let a = vec![vec![2, 1, 1], vec![3, 1, 2], vec![2, 1, 2]];

        // Using your original B values:
        let b = vec![8, 11, 3];

        let result = gaussian_elimination(a, b);

        // The actual mathematical solution for this specific system is:
        let expected = vec![8.0, -3.0, -5.0];

        // Use a tolerance check instead of assert_eq!
        let epsilon = 1e-10;
        for i in 0..result.len() {
            assert!(
                (result[i] - expected[i]).abs() < epsilon,
                "Values at index {} do not match: got {}, expected {}",
                i,
                result[i],
                expected[i]
            );
        }
    }
}