Modified to upmove the row with highest pivot & added more comments

Author: rishvic

contents/gaussian_elimination/code/c++/gaussian_elimination.cpp

@@ -2,6 +2,7 @@
 #include <algorithm>
 #include <stdexcept>
 #include <vector>
+#include <cmath>
 using namespace std;
 
 class EquationSolver
@@ -29,41 +30,45 @@ class EquationSolver
   // to swap.
 
 
-  void gaussianElimination(vector<vector<long double>> &eqns, int varc)
+  int gaussianElimination(vector<vector<long double>> &eqns, int varc)
   {
     // 'eqns' is the matrix, 'varc' is no. of vars
+    int err = 0; // marker for if matrix is singular
     for (int i = 0; i < varc - 1; i++)
     {
-      if (eqns[i][i] == 0)
+      long double pivot = fabsl(eqns[i][i]);
+      int newRow = i;
+      for (int j = i + 1; j < varc; j++)
       {
-        for (int j = i + 1; j < varc; j++)
+        if (fabsl(eqns[j][i]) > pivot)
         {
-          if (eqns[j][i] != 0)
-          {
-            swapRow(eqns[j], eqns[i], i);
-            break;
-          }
+          pivot = fabsl(eqns[j][i]);
+          newRow = j;
         }
+      }
 
-        if (eqns[i][i] == 0)
-          throw runtime_error("Singular matrix");
-        // Cannot solve since coefficient of one variable is always zero
+      if (pivot == 0.0)
+      {
+        err = 1; // Marking that matrix is singular
+        continue;
       }
+      if (newRow != i)
+        swapRow(eqns[newRow], eqns[i], i);
 
       for (int j = i + 1; j < varc; j++)
         if (eqns[j][i])
           subRow(eqns[j], eqns[i], (eqns[j][i] / eqns[i][i]), i);
     }
 
-    if (eqns[varc - 1][varc - 1] == 0)
+    if (eqns[varc - 1][varc - 1] == 0 || err)
     {
-      if (eqns[varc - 1][varc] == 0)
-        throw runtime_error("Singular matrix");
-      // Cannot solve since final equation is '0*xn = 0'
+      if (eqns[varc - 1][varc] == 0 || err)
+        return 1; // Error code: Singular matrix
       else
-        throw runtime_error("No solutions");
+        return 2; // Error code: No solutions
       // Cannot solve since final equation is '0*xn = c'
     }
+    return 0; // Successful
   }
 
 
@@ -73,11 +78,12 @@ class EquationSolver
     for (int i = varc - 1; i >= 0; i--)
     {
       eqns[i][varc] /= eqns[i][i];
-      eqns[i][i] = 1;
-      for (int j = i - 1; j >= 0; j--)
+      eqns[i][i] = 1; // We know that the only entry in this row is 1
+
+      for (int j = i - 1; j >= 0; j--) // subtracting rows from below
       {
         eqns[j][varc] -= eqns[j][i] * eqns[i][varc];
-        eqns[j][i] = 0;
+        eqns[j][i] = 0; // We also set all the other values in row to 0 directly
       }
     }
   }
@@ -101,7 +107,19 @@ class EquationSolver
 public:
   vector<long double> solveByGaussJordan(vector< vector<long double> > eqns, int varc)
   {
-    this->gaussianElimination(eqns, varc);
+    int status = this->gaussianElimination(eqns, varc);
+    switch (status)
+    {
+      case 0:
+        break;
+
+      case 1:
+        throw runtime_error("Singular matrix");
+
+      case 2:
+        throw runtime_error("No solutions");
+    }
+
     this->gaussJordan(eqns, varc);
 
     vector<long double> ans(varc);
@@ -112,30 +130,57 @@ class EquationSolver
 
   vector<long double> solveByBacksubs(vector< vector<long double> > eqns, int varc)
   {
-    this->gaussianElimination(eqns, varc);
+    int status = this->gaussianElimination(eqns, varc);
+    switch (status)
+    {
+      case 0:
+        break;
+
+      case 1:
+        throw runtime_error("Singular matrix");
+
+      case 2:
+        throw runtime_error("No solutions");
+    }
+
     vector<long double> ans = this->backSubs(eqns, varc);
     return ans;
   }
 };
 
 int main() {
   int varc = 3;
-  vector< vector<long double> > equations { { 2 , 3, 4, 6 },
-                                            { 1 , 2, 3, 4 },
-                                            { 3 ,-4, 0, 10 }};
+  vector< vector<long double> > equations { { 2, 3, 4, 6 },
+                                            { 1, 2, 3, 4 },
+                                            { 3, -4, 0, 10 }};
   EquationSolver solver;
   vector<long double> ans;
   try
   {
     ans = solver.solveByGaussJordan(equations, varc);
   }
-  catch(runtime_error e)
+  catch(runtime_error &e)
+  {
+    cout << "Error found: " << e.what() << endl;
+    return -1;
+  }
+
+  cout << "The solution is (by Gauss-Jordan)," << endl
+       << "x == " << ans[0] << endl
+       << "y == " << ans[1] << endl
+       << "z == " << ans[2] << endl << endl;
+
+  try
+  {
+    ans = solver.solveByBacksubs(equations, varc);
+  }
+  catch (runtime_error &e)
   {
     cout << "Error found: " << e.what() << endl;
     return -1;
   }
 
-  cout << "The solution is," << endl
+  cout << "The solution is (by Backsubstitution)," << endl
        << "x == " << ans[0] << endl
        << "y == " << ans[1] << endl
        << "z == " << ans[2] << endl;

contents/gaussian_elimination/gaussian_elimination.md

@@ -315,8 +315,8 @@ In code, this process might look like this:
 [import:5-13, lang:"c_cpp"](code/c/gaussian_elimination.c)
 [import:19-34, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="cpp" %}
-[import:22-29, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
-[import:37-51, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
+[import:23-30, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
+[import:36-56, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
 {% sample lang="hs" %}
 [import:10-17, lang:"haskell"](code/haskell/gaussianElimination.hs)
 [import:44-46, lang:"haskell"](code/haskell/gaussianElimination.hs)
@@ -388,7 +388,7 @@ Here is what it might look like in code:
 {% sample lang="c" %}
 [import:36-41, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="cpp" %}
-[import:11-19, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
+[import:12-20, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
 {% sample lang="hs" %}
 [import:19-33, lang:"haskell"](code/haskell/gaussianElimination.hs)
 [import:42-42, lang:"haskell"](code/haskell/gaussianElimination.hs)
@@ -409,7 +409,7 @@ When we put everything together, it looks like this:
 {% sample lang="c" %}
 [import:15-48, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="cpp" %}
-[import:11-67, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
+[import:12-72, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
 {% sample lang="rs" %}
 [import:41-78, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% sample lang="hs" %}
@@ -448,7 +448,7 @@ Here it is in code:
 {% sample lang="c" %}
 [import:64-82, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="cpp" %}
-[import:70-83, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
+[import:75-89, lang:"c_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)
@@ -491,7 +491,7 @@ In code, it looks like this:
 {% sample lang="c" %}
 [import:50-62, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="cpp" %}
-[import:86-99, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
+[import:92-105, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
 {% sample lang="rs" %}
 [import:79-94, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% sample lang="hs" %}