Removed singular matrix checking & made the code simpler

Author: rishvic

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

@@ -5,60 +5,49 @@
 #include <vector>
 
 
-int gaussianElimination(std::vector< std::vector<double> > &eqns, int varc) {
+void gaussianElimination(std::vector< std::vector<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++) {
-    double pivot = fabs(eqns[i][i]);
-    int newRow = i;
+    int pivot = i;
 
-    for (int j = i + 1; j < varc; j++) {
-      if (fabs(eqns[j][i]) > pivot) {
-        pivot = fabs(eqns[j][i]);
-        newRow = j;
-      }
-    }
+    for (int j = i + 1; j < varc; j++)
+      if (fabs(eqns[j][i]) > fabs(eqns[pivot][i]))
+        pivot = j;
 
-    if (pivot == 0.0) {
+    if (eqns[pivot][i] == 0.0) {
       err = 1; // Marking that matrix is singular
-      continue;
+      continue; // But continuing to simplify the matrix as much as possible
     }
 
-    if (newRow != i) // Swapping the rows if new row with higher pivot is found
-      std::swap(eqns[newRow], eqns[i]); // C++ swap function
+    if (i != pivot) // Swapping the rows if new row with higher maxVals is found
+      std::swap(eqns[pivot], eqns[i]); // C++ swap function
 
     for (int j = i + 1; j < varc; j++) {
-      if (eqns[j][i]) {
-        double mult = eqns[j][i] / eqns[i][i];
-        eqns[j][i] = 0.0;
-        for (int k = i + 1; k < (varc + 1); k++) // k doesn't start at 0, since
-          eqns[j][k] -= mult * eqns[i][k];       // values before from 0 to i
-      }                                          // are already 0
-    }
-  }
+      double scale = eqns[j][i] / eqns[i][i];
 
-  if (eqns[varc - 1][varc - 1] == 0 || err) {
-    if (eqns[varc - 1][varc] == 0 || err)
-      return 1; // Error code: Singular matrix
-    else
-      return 2; // Error code: No solutions
-    // Cannot solve since final equation is '0*xn = c'
+      for (int k = i + 1; k < (varc + 1); k++) // k doesn't start at 0, since
+        eqns[j][k] -= scale * eqns[i][k];      // values before from 0 to i
+                                               // are already 0
+      eqns[j][i] = 0.0;
+    }
   }
-
-  return 0; // Successful
 }
 
 
 void gaussJordan(std::vector< std::vector<double> > &eqns, int varc) {
   // 'eqns' is the (Row-echelon) matrix, 'varc' is no. of vars
   for (int i = varc - 1; i >= 0; i--) {
-    eqns[i][varc] /= eqns[i][i];
-    eqns[i][i] = 1; // We know that the only entry in this row is 1
-
-    // subtracting rows from below
-    for (int j = i - 1; j >= 0; j--) {
-      eqns[j][varc] -= eqns[j][i] * eqns[i][varc];
-      eqns[j][i] = 0; // We also set all the other values in row to 0 directly
+    if (eqns[i][i] != 0) {
+      eqns[i][varc] /= eqns[i][i];
+      eqns[i][i] = 1; // We know that the only entry in this row is 1
+
+      // subtracting rows from below
+      for (int j = i - 1; j >= 0; j--) {
+        eqns[j][varc] -= eqns[j][i] * eqns[i][varc];
+        eqns[j][i] = 0; // We also set all the other values in row to 0 directly
+      }
     }
   }
 }
@@ -69,30 +58,22 @@ std::vector<double> backSubs(std::vector<std::vector<double>> &eqns, int varc)
   // 'eqns' is matrix, 'varc' is no. of variables
   std::vector<double> ans(varc);
   for (int i = varc - 1; i >= 0; i--) {
-    double sum = 0;
+    double sum = 0.0;
+
     for (int j = i + 1; j < varc; j++)
       sum += eqns[i][j] * ans[j];
 
-    ans[i] = (eqns[i][varc] - sum) / eqns[i][i];
+    if (eqns[i][i] != 0)
+      ans[i] = (eqns[i][varc] - sum) / eqns[i][i];
+    else
+      return std::vector<double>(0);
   }
   return ans;
 }
 
 
 std::vector<double> solveByGaussJordan(std::vector<std::vector<double>> eqns, int varc) {
-  int status = gaussianElimination(eqns, varc);
-  switch (status) {
-  case 0:
-    break;
-
-  case 1:
-    std::cout << "Singular matrix" << std::endl;
-    return std::vector<double>(0); // return empty vector in case of no solutions
-
-  case 2:
-    std::cout << "No solutions" << std::endl;
-    return std::vector<double>(0); // return empty vectors in case of no solutions
-  }
+  gaussianElimination(eqns, varc);
 
   gaussJordan(eqns, varc);
 
@@ -104,19 +85,7 @@ std::vector<double> solveByGaussJordan(std::vector<std::vector<double>> eqns, in
 
 
 std::vector<double> solveByBacksubs(std::vector<std::vector<double>> eqns, int varc) {
-  int status = gaussianElimination(eqns, varc);
-  switch (status) {
-  case 0:
-    break;
-
-  case 1:
-    std::cout << "Singular matrix" << std::endl;
-    return std::vector<double>(0); // return empty vectors in case of no solutions
-
-  case 2:
-    std::cout << "No solutions" << std::endl;
-    return std::vector<double>(0); // return empty vectors in case of no solutions
-  }
+  gaussianElimination(eqns, varc);
 
   std::vector<double> ans = backSubs(eqns, varc);
   return ans;
@@ -134,23 +103,15 @@ int main() {
 
   ans = solveByGaussJordan(equations, varc);
 
-  if (ans.empty()) { // If ans is empty, then no unique solutions exist
-    std::cout << "No solution" << std::endl;
-  } else {
-    std::cout << "The solution is (by Gauss Jordan)," << std::endl
-         << "x == " << ans[0] << std::endl
-         << "y == " << ans[1] << std::endl
-         << "z == " << ans[2] << std::endl;
-  }
+  std::cout << "The solution is (by Gauss Jordan)," << std::endl
+        << "x == " << ans[0] << std::endl
+        << "y == " << ans[1] << std::endl
+        << "z == " << ans[2] << std::endl;
 
   ans = solveByBacksubs(equations, varc);
 
-  if (ans.empty()) { // If ans is empty, then no unique solutions exist
-    std::cout << "No solution" << std::endl;
-  } else {
-    std::cout << "The solution is (by Backsubstitution)," << std::endl
-         << "x == " << ans[0] << std::endl
-         << "y == " << ans[1] << std::endl
-         << "z == " << ans[2] << std::endl;
-  }
+  std::cout << "The solution is (by Backsubstitution)," << std::endl
+        << "x == " << ans[0] << std::endl
+        << "y == " << ans[1] << std::endl
+        << "z == " << ans[2] << std::endl;
 }

contents/gaussian_elimination/gaussian_elimination.md

@@ -315,8 +315,8 @@ In code, this process might look like this:
 [import:5-13, lang:"c"](code/c/gaussian_elimination.c)
 [import:19-34, lang:"c"](code/c/gaussian_elimination.c)
 {% sample lang="cpp" %}
-[import:27-28, lang:"cpp"](code/c++/gaussian_elimination.cpp)
-[import:12-28, lang:"cpp"](code/c++/gaussian_elimination.cpp)
+[import:25-25, lang:"cpp"](code/c++/gaussian_elimination.cpp)
+[import:13-25, lang:"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"](code/c/gaussian_elimination.c)
 {% sample lang="cpp" %}
-[import:32-35, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
+[import:28-32, 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"](code/c/gaussian_elimination.c)
 {% sample lang="cpp" %}
-[import:8-49, lang:"cpp"](code/c++/gaussian_elimination.cpp)
+[import:8-36, lang:"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"](code/c/gaussian_elimination.c)
 {% sample lang="cpp" %}
-[import:52-64, lang:"cpp"](code/c++/gaussian_elimination.cpp)
+[import:39-53, 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)
@@ -491,7 +491,7 @@ In code, it looks like this:
 {% sample lang="c" %}
 [import:50-62, lang:"c"](code/c/gaussian_elimination.c)
 {% sample lang="cpp" %}
-[import:67-79, lang:"cpp"](code/c++/gaussian_elimination.cpp)
+[import:56-72, lang:"cpp"](code/c++/gaussian_elimination.cpp)
 {% sample lang="rs" %}
 [import:79-94, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% sample lang="hs" %}