Implement Gaussian Elimination in C++

Author: mika314

chapters/algorithms/gaussian_elimination/code/c++/gaussian_elimination.cpp

@@ -0,0 +1,99 @@
+#include <algorithm>
+#include <cassert>
+#include <iostream>
+#include <vector>
+
+void gaussian_elimination(std::vector<double>& a, size_t cols) {
+  assert(a.size() % cols == 0);
+  auto rows = a.size() / cols;
+
+  // Main loop going through all columns
+  for (auto k = 0u; k < std::min(rows, cols); ++k) {
+    // Step 1: finding the maximum element for each column
+    auto max_index = [&]() {
+      auto res = k;
+      auto max_element = a[k + k * cols];
+      for (auto r = k + 1; r < rows; ++r)
+        if (max_element < a[k + r * cols]) {
+          max_element = a[k + r * cols];
+          res = r;
+        }
+      return res;
+    }();
+
+    // Check to make sure matrix is good!
+    if (a[k + max_index * cols] == 0) {
+      std::cout << "matrix is singular! End!\n";
+      return;
+    }
+
+    // Step 2: swap row with highest value for that column to the top
+    for (auto c = 0u; c < cols; ++c)
+      std::swap(a[c + k * cols], a[c + max_index * cols]);
+
+    // Loop for all remaining rows
+    for (auto i = k + 1; i < rows; ++i) {
+
+      // Step 3: finding fraction
+      auto fraction = a[k + i * cols] / a[k + k * cols];
+
+      // loop through all columns for that row
+      for (auto j = k + 1; j < cols; ++j) {
+
+        // Step 4: re-evaluate each element
+        a[j + i * cols] -= a[j + k * cols] * fraction;
+      }
+
+      // Step 5: Set lower elements to 0
+      a[k + i * cols] = 0;
+    }
+  }
+}
+
+std::vector<double> back_substitution(std::vector<double>& a, size_t cols) {
+  assert(a.size() % cols == 0);
+  auto rows = a.size() / cols;
+
+  // Creating the solution Vector
+  std::vector<double> soln;
+  soln.resize(rows);
+
+  // initialize the final element
+  soln[rows - 1] =
+      a[cols - 1 + (rows - 1) * cols] / a[cols - 1 - 1 + (rows - 1) * cols];
+
+  for (int i = rows - 1; i >= 0; --i) {
+    auto sum = 0.0;
+    for (int j = cols - 2; j > i; --j) {
+      sum += soln[j] * a[j + i * cols];
+    }
+
+    soln[i] = (a[cols - 1 + i * cols] - sum) / a[i + i * cols];
+  }
+
+  return soln;
+}
+
+int main() {
+  std::vector<double> a = {2,  3,  4,  6,
+                           1,  2,  3,  4,
+                           3, -4,  0, 10};
+  const auto cols = 4u;
+  assert(a.size() % cols == 0);
+  const auto rows = a.size() / cols;
+
+  gaussian_elimination(a, cols);
+
+  for (auto i = 0u; i < rows; ++i) {
+    std::cout << "[";
+    for (auto j = 0u; j < cols; ++j) {
+      std::cout << a[j + i * cols] << " ";
+    }
+    std::cout << "]\n";
+  }
+
+  auto soln = back_substitution(a, 4);
+
+  for (auto element : soln)
+    std::cout << element << std::endl;
+}

chapters/algorithms/gaussian_elimination/gaussian_elimination.md

@@ -259,6 +259,8 @@ In code, this looks like:
 [import:1-42, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
 [import:13-44, lang:"c_cpp"](code/c/gaussian_elimination.c)
+{% sample lang="cpp" %}
+[import:6-51, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
 {% sample lang="rs" %}
 [import:41-78, lang:"rust"](code/rust/gaussian_elimination.rs)
 {% endmethod %}