"Gaussian Elimination" Implementation in C++

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

@@ -0,0 +1,114 @@
+#include <algorithm>
+#include <cmath>
+#include <iomanip>
+#include <iostream>
+#include <vector>
+
+
+void gaussianElimination(std::vector<std::vector<double> > &eqns) {
+  // 'eqns' is the matrix, 'rows' is no. of vars
+  int rows = eqns.size(), cols = eqns[0].size();
+
+  for (int i = 0; i < rows - 1; i++) {
+    int pivot = i;
+
+    for (int j = i + 1; j < rows; j++) {
+      if (fabs(eqns[j][i]) > fabs(eqns[pivot][i])) pivot = j;
+    }
+
+    if (eqns[pivot][i] == 0.0)
+      continue;  // But continuing to simplify the matrix as much as possible
+
+    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 < rows; j++) {
+      double scale = eqns[j][i] / eqns[i][i];
+
+      for (int k = i + 1; k < cols; 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;
+    }
+  }
+}
+
+void gaussJordan(std::vector<std::vector<double> > &eqns) {
+  // 'eqns' is the (Row-echelon) matrix, 'rows' is no. of vars
+  int rows = eqns.size();
+
+  for (int i = rows - 1; i >= 0; i--) {
+
+    if (eqns[i][i] != 0) {
+
+      eqns[i][rows] /= 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][rows] -= eqns[j][i] * eqns[i][rows];
+        eqns[j][i] = 0;  // We also set all the other values in row to 0 directly
+      }
+    }
+  }
+}
+
+std::vector<double> backSubs(const std::vector<std::vector<double> > &eqns) {
+  // 'eqns' is matrix, 'rows' is no. of variables
+  int rows = eqns.size();
+
+  std::vector<double> ans(rows);
+  for (int i = rows - 1; i >= 0; i--) {
+    double sum = 0.0;
+
+    for (int j = i + 1; j < rows; j++) sum += eqns[i][j] * ans[j];
+
+    if (eqns[i][i] != 0)
+      ans[i] = (eqns[i][rows] - sum) / eqns[i][i];
+    else
+      return std::vector<double>(0);
+  }
+  return ans;
+}
+
+
+void printMatrix(const std::vector<std::vector<double> > &matrix) {
+  for (int row = 0; row < matrix.size(); row++) {
+    std::cout << "[";
+
+    for (int col = 0; col < matrix[row].size() - 1; col++)
+      std::cout << std::setw(8) << std::fixed << std::setprecision(3)
+                << matrix[row][col];
+
+    std::cout << " |" << std::setw(8) << std::fixed << std::setprecision(3)
+              << matrix[row].back() << " ]" << std::endl;
+  }
+}
+
+
+int main() {
+  std::vector<std::vector<double> > equations{
+      {2, 3, 4, 6},
+      {1, 2, 3, 4},
+      {3, -4, 0, 10}};
+
+  std::cout << "Initial matrix:" << std::endl;
+  printMatrix(equations);
+  std::cout << std::endl;
+
+  gaussianElimination(equations);
+  std::cout << "Matrix after gaussian elimination:" << std::endl;
+  printMatrix(equations);
+  std::cout << std::endl;
+
+  std::vector<double> ans = backSubs(equations);
+  std::cout << "Solution from backsubstitution" << std::endl;
+  std::cout << "x = " << ans[0] << ", y = " << ans[1] << ", z = " << ans[2]
+            << std::endl
+            << std::endl;
+
+  gaussJordan(equations);
+  std::cout << "Matrix after Gauss Jordan:" << std::endl;
+  printMatrix(equations);
+  std::cout << std::endl;
+}

contents/gaussian_elimination/gaussian_elimination.md

@@ -314,6 +314,8 @@ In code, this process might look like this:
 {% sample lang="c" %}
 [import:5-13, lang:"c"](code/c/gaussian_elimination.c)
 [import:19-34, lang:"c"](code/c/gaussian_elimination.c)
+{% sample lang="cpp" %}
+[import:13-23, 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)
@@ -384,6 +386,8 @@ Here is what it might look like in code:
 [import:32-40, lang:"java"](code/java/GaussianElimination.java)
 {% sample lang="c" %}
 [import:36-41, lang:"c"](code/c/gaussian_elimination.c)
+{% sample lang="cpp" %}
+[import:25-32, lang:"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)
@@ -403,6 +407,8 @@ When we put everything together, it looks like this:
 [import:1-45, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
 [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" %}
@@ -440,6 +446,8 @@ Here it is in code:
 [import:67-93, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
 [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)
@@ -481,6 +489,8 @@ In code, it looks like this:
 [import:47-64, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
 [import:50-62, lang:"c"](code/c/gaussian_elimination.c)
+{% 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)
 {% sample lang="hs" %}
@@ -547,6 +557,8 @@ Here's a video describing Gaussian elimination:
 [import, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
 [import, lang:"c"](code/c/gaussian_elimination.c)
+{% sample lang="cpp" %}
+[import, lang:"cpp"](code/c++/gaussian_elimination.cpp)
 {% sample lang="rs" %}
 [import, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% sample lang="hs" %}