Gauss Jordan elimination in rust (algorithm-archivists#559)

Author: Gathros

contents/gaussian_elimination/code/rust/gaussian_elimination.rs

@@ -76,6 +76,25 @@ fn gaussian_elimination(a: &mut Matrix) {
     }
 }
 
+fn gauss_jordan(a: &mut Matrix) {
+    let mut row = 0;
+    for k in 0..(a.cols - 1) {
+        if a[(row, k)] != 0.0 {
+            for i in (k..a.cols).rev() {
+                a[(row, i)] /= a[(row, k)];
+            }
+
+            for i in 0..row {
+                for j in (k..a.cols).rev() {
+                    a[(i, j)] -= a[(i, k)] * a[(row, j)];
+                }
+            }
+
+            row += 1;
+        }
+    }
+}
+
 fn back_substitution(a: &Matrix) -> Vec<f64> {
     let mut soln = vec![0.0; a.rows];
 
@@ -101,6 +120,7 @@ fn main() {
     );
 
     gaussian_elimination(&mut a);
+    gauss_jordan(&mut a);
     let soln = back_substitution(&a);
     println!("Solution: {:?}", soln);
 }

contents/gaussian_elimination/gaussian_elimination.md

@@ -324,7 +324,7 @@ In code, this process might look like this:
 {% sample lang="py" %}
 [import:13-19, lang:"python"](code/python/gaussian_elimination.py)
 {% sample lang="rs" %}
-[import:43-60, lang:"rust"](code/rust/gaussian_elimination.rs)
+[import:43-76, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% endmethod %}
 
 As a note, if the highest value is $$0$$, the matrix is singular and the system has no single solution.
@@ -396,7 +396,7 @@ Here is what it might look like in code:
 {% sample lang="py" %}
 [import:21-26, lang:"python"](code/python/gaussian_elimination.py)
 {% sample lang="rs" %}
-[import:62-71, lang:"rust"](code/rust/gaussian_elimination.rs)
+[import:62-75, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% endmethod %}
 
 #### All together
@@ -409,8 +409,6 @@ When we put everything together, it looks like this:
 [import:15-48, lang:"c"](code/c/gaussian_elimination.c)
 {% sample lang="cpp" %}
 [import:8-34, lang:"cpp"](code/c++/gaussian_elimination.cpp)
-{% sample lang="rs" %}
-[import:41-78, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% sample lang="hs" %}
 [import:10-36, lang:"haskell"](code/haskell/gaussianElimination.hs)
 {% sample lang="py" %}
@@ -421,6 +419,8 @@ When we put everything together, it looks like this:
 [import:1-38, lang:"javascript"](code/javascript/gaussian_elimination.js)
 {% sample lang="go" %}
 [import:9-53, lang:"go"](code/go/gaussian_elimination.go)
+{% sample lang="rs" %}
+[import:41-77, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% endmethod %}
 
 To be clear: if the matrix is found to be singular during this process, the system of equations is either over- or under-determined and no general solution exists.
@@ -450,9 +450,6 @@ Here it is in code:
 [import:64-82, lang:"c"](code/c/gaussian_elimination.c)
 {% sample lang="cpp" %}
 [import:36-54, lang:"cpp"](code/c++/gaussian_elimination.cpp)
-{% sample lang="rs" %}
-This code does not exist yet in rust, so here's Julia code (sorry for the inconvenience)
-[import:67-93, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="hs" %}
 [import:38-46, lang:"haskell"](code/haskell/gaussianElimination.hs)
 {% sample lang="py" %}
@@ -463,6 +460,8 @@ This code does not exist yet in rust, so here's Julia code (sorry for the inconv
 [import:57-76, lang:"javascript"](code/javascript/gaussian_elimination.js)
 {% sample lang="go" %}
 [import:55-82, lang:"go"](code/go/gaussian_elimination.go)
+{% sample lang="rs" %}
+[import:79-96, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% endmethod %}
 
 As a note: Gauss-Jordan elimination can also be used to find the inverse of a matrix by following the same procedure to generate a reduced row echelon matrix, but with an identity matrix on the other side instead of the right-hand side of each equation.
@@ -496,7 +495,7 @@ In code, it looks like this:
 {% sample lang="cpp" %}
 [import:56-72, lang:"cpp"](code/c++/gaussian_elimination.cpp)
 {% sample lang="rs" %}
-[import:79-94, lang:"rust"](code/rust/gaussian_elimination.rs)
+[import:98-112, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% sample lang="hs" %}
 [import:48-53, lang:"haskell"](code/haskell/gaussianElimination.hs)
 {% sample lang="py" %}
@@ -523,7 +522,7 @@ That said, it is particularly interesting to see what happens as we plot our mat
 <video width="560" height="315" autoplay controls loop>
   <source src="res/GE_vis.mp4" type="video/mp4">
 Your browser does not support the video tag.
-</video> 
+</video>
 </div>
 
 As we can see in the above visualization, the planes wobble about in 3D until they reach row echelon form, where one plane is parallel to the $$x$$ and $$y$$ axes.