Skip to main content

solver/linalg/
gaussian.rs

1/// Simple Gaussian Elimination solver for dense linear systems.
2/// Solves $Ax = b$ where $A$ is an $n \times n$ matrix provided as a flat row-major slice.
3///
4/// # Arguments
5/// * `a` - The coefficient matrix $A$ (flat slice of size $n^2$).
6/// * `b` - The right-hand side vector $b$ (slice of size $n$).
7/// * `n` - The dimension of the system.
8///
9/// # Returns
10/// A `Vec<f64>` containing the solution $x$.
11pub fn solve_linear_system(a: &[f64], b: &[f64], n: usize) -> Vec<f64> {
12    let mut aug = vec![0.0; n * (n + 1)];
13    // Flatten A|b into augmented matrix
14    for i in 0..n {
15        for j in 0..n {
16            aug[i * (n + 1) + j] = a[i * n + j];
17        }
18        aug[i * (n + 1) + n] = b[i];
19    }
20
21    // Forward elimination with pivoting
22    for i in 0..n {
23        let mut pivot = i;
24        for j in i + 1..n {
25            if aug[j * (n + 1) + i].abs() > aug[pivot * (n + 1) + i].abs() {
26                pivot = j;
27            }
28        }
29
30        if pivot != i {
31            // Swap rows
32            for k in 0..=n {
33                // =n because it's augmented
34                aug.swap(i * (n + 1) + k, pivot * (n + 1) + k);
35            }
36        }
37
38        if aug[i * (n + 1) + i].abs() < 1e-12 {
39            continue; // Singular or near-singular
40        }
41
42        for j in i + 1..n {
43            let factor = aug[j * (n + 1) + i] / aug[i * (n + 1) + i];
44            for k in i..=n {
45                aug[j * (n + 1) + k] -= factor * aug[i * (n + 1) + k];
46            }
47        }
48    }
49
50    // Back substitution
51    let mut x = vec![0.0; n];
52    for i in (0..n).rev() {
53        let mut sum = 0.0;
54        for j in i + 1..n {
55            sum += aug[i * (n + 1) + j] * x[j];
56        }
57        if aug[i * (n + 1) + i].abs() > 1e-12 {
58            x[i] = (aug[i * (n + 1) + n] - sum) / aug[i * (n + 1) + i];
59        }
60    }
61    x
62}