solver/linalg/
gaussian.rs1pub fn solve_linear_system(a: &[f64], b: &[f64], n: usize) -> Vec<f64> {
12 let mut aug = vec![0.0; n * (n + 1)];
13 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 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 for k in 0..=n {
33 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; }
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 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}