"Gaussian Elimination" Implementation in C++

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

@@ -0,0 +1,117 @@
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include <stdexcept>
+#include <vector>
+
+
+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++) {
+ int pivot = i;
+
+ for (int j = i + 1; j < varc; j++)
+ if (fabs(eqns[j][i]) > fabs(eqns[pivot][i]))
+ pivot = j;
+
+ if (eqns[pivot][i] == 0.0) {
+ err = 1; // Marking that matrix is singular
+ 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 < varc; j++) {
+ double scale = eqns[j][i] / eqns[i][i];
+
+ 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;
+ }
+ }
+}
+
+
+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--) {
+ 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
+ }
+ }
+ }
+}
+
+
+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.0;
+
+ for (int j = i + 1; j < varc; j++)
+ sum += eqns[i][j] * ans[j];
+
+ 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) {
+ gaussianElimination(eqns, varc);
+
+ gaussJordan(eqns, varc);
+
+ std::vector<double> ans(varc);
+ for (int i = 0; i < varc; i++)
+ ans[i] = eqns[i][varc];
+ return ans;
+}
+
+
+std::vector<double> solveByBacksubs(std::vector<std::vector<double>> eqns, int varc) {
+ gaussianElimination(eqns, varc);
+
+ std::vector<double> ans = backSubs(eqns, varc);
+ return ans;
+}
+
+
+int main() {
+ int varc = 3;
+ std::vector< std::vector<double> > equations{
+ {2, 3, 4, 6},
+ {1, 2, 3, 4},
+ {3, -4, 0, 10}};
+
+ std::vector<double> ans;
+
+ ans = solveByGaussJordan(equations, varc);
+
+ 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);
+
+ 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

@@ -232,7 +232,7 @@ $$
 $$
 
 There are plenty of different strategies you could use to do this, and no one strategy is better than the rest.
-One method is to subtract a multiple of the top row from subsequent rows below it such that all values beneath the pivot value are zero. 
+One method is to subtract a multiple of the top row from subsequent rows below it such that all values beneath the pivot value are zero.
 This process might be easier if you swap some rows around first and can be performed for each pivot.
 
 After you get a row echelon matrix, the next step is to find the reduced row echelon form. In other words, we do the following:
@@ -312,8 +312,11 @@ In code, this process might look like this:
 {% sample lang="java" %}
 [import:14-30, lang:"java"](code/java/GaussianElimination.java)
 {% sample lang="c" %}
-[import:5-13, lang:"c_cpp"](code/c/gaussian_elimination.c)
-[import:19-34, lang:"c_cpp"](code/c/gaussian_elimination.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: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)
@@ -367,7 +370,7 @@ $$
 \left[
 \begin{array}{ccc|c}
 3 & -4 & 0 & 10 \\
-0 & \mathbf{\frac{10}{3}} & \mathbf{3} & \mathbf{\frac{2}{3}} 
+0 & \mathbf{\frac{10}{3}} & \mathbf{3} & \mathbf{\frac{2}{3}}
 \\
 2 & 3 & 4 & 6
 \end{array}
@@ -383,7 +386,9 @@ Here is what it might look like in code:
 {% sample lang="java" %}
 [import:32-40, lang:"java"](code/java/GaussianElimination.java)
 {% sample lang="c" %}
-[import:36-41, lang:"c_cpp"](code/c/gaussian_elimination.c)
+[import:36-41, lang:"c"](code/c/gaussian_elimination.c)
+{% sample lang="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)
@@ -402,7 +407,9 @@ When we put everything together, it looks like this:
 {% sample lang="jl" %}
 [import:1-45, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
-[import:15-48, lang:"c_cpp"](code/c/gaussian_elimination.c)
+[import:15-48, lang:"c"](code/c/gaussian_elimination.c)
+{% sample lang="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" %}
@@ -439,7 +446,9 @@ Here it is in code:
 {% sample lang="jl" %}
 [import:67-93, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
-[import:64-82, lang:"c_cpp"](code/c/gaussian_elimination.c)
+[import:64-82, lang:"c"](code/c/gaussian_elimination.c)
+{% sample lang="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)
@@ -480,7 +489,9 @@ In code, it looks like this:
 {% sample lang="jl" %}
 [import:47-64, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
-[import:50-62, lang:"c_cpp"](code/c/gaussian_elimination.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" %}
@@ -495,7 +506,7 @@ In code, it looks like this:
 
 ## Visual Representation
 
-We have thus far used Gaussian elimination as a method to solve a system of equations; however, there is often a much easier way to find a similar solution simply by plotting each row in our matrix. 
+We have thus far used Gaussian elimination as a method to solve a system of equations; however, there is often a much easier way to find a similar solution simply by plotting each row in our matrix.
 For the case of 2 equations and 2 unknowns, we would plot the two lines corresponding to each equation and the $$(x, y)$$ location of their point of intersection would be the solution for $$x$$ and $$y$$.
 Similarly, for the case of 3 equations and 3 unknowns, we would plot 3 planes and the $$(x, y, z)$$ location of their point of intersection would be the solution for $$x$$, $$y$$, and $$z$$.
 
@@ -535,7 +546,9 @@ There are also plenty of other solvers that do similar things that we will get t
 {% sample lang="jl" %}
 [import, lang:"julia"](code/julia/gaussian_elimination.jl)
 {% sample lang="c" %}
-[import, lang:"c_cpp"](code/c/gaussian_elimination.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" %}