Skip to content

"Gaussian Elimination" Implementation in C++ #555

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 9 commits into from
Jan 7, 2019
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
187 changes: 187 additions & 0 deletions contents/gaussian_elimination/code/c++/gaussian_elimination.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,187 @@
#include <iostream>
#include <algorithm>
#include <stdexcept>
#include <vector>
#include <cmath>
using namespace std;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think the use of using namespace std is a good idea for the philosophy of AAA is. The idea is for the AAA to be used by new comers to the language to learn, also for people to use the algorithm in real work cases and show case the algorithm. I don't believe that using namespace std achieves this.


class EquationSolver
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There really is no need in making a class for solving Gaussian elimination.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So, should I just define each function individually, because I thought encapsulating it would make it more user-friendly and show how it can be used to solve equations.

{
private:

// subtract row from another row after multiplcation
void subRow(vector<long double> &toBeSubtracted, const vector<long double> &toSubtract,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is no need of making this function if it's used in one place.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since the row subtraction function was being showcased in the documentation, I thought a special mention would be necessary.

long double mult, int start)
{
for (int i = start; i < toBeSubtracted.size(); i++)
toBeSubtracted[i] -= mult * toSubtract[i];
}
// A start variable is given since if we know all values are 0, we don't need
// to subtract.


// swap two rows
void swapRow(vector<long double> &eqn1, vector<long double> &eqn2, int start)
{
for (int i = start; i < eqn1.size(); i++)
swap(eqn1[i], eqn2[i]);
}
// A start variable is given since if we know all values are 0, we don't need
// to swap.


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++)
{
long double pivot = fabsl(eqns[i][i]);
int newRow = i;
for (int j = i + 1; j < varc; j++)
{
if (fabsl(eqns[j][i]) > pivot)
{
pivot = fabsl(eqns[j][i]);
newRow = j;
}
}

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 || 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'
}
return 0; // Successful
}


void gaussJordan(vector<vector<long 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

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; // We also set all the other values in row to 0 directly
}
}
}


vector<long double> backSubs(vector<vector<long double>> &eqns, int varc)
{
// 'eqns' is matrix, 'varc' is no. of variables
vector<long double> ans(varc);
for (int i = varc - 1; i >= 0; i--)
{
long double sum = 0;
for (int j = i + 1; j < varc; j++)
sum += eqns[i][j] * ans[j];

ans[i] = (eqns[i][varc] - sum) / eqns[i][i];
}
return ans;
}

public:
vector<long double> solveByGaussJordan(vector< vector<long double> > eqns, int 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);
for (int i = 0; i < varc; i++)
ans[i] = eqns[i][varc];
return ans;
}

vector<long double> solveByBacksubs(vector< vector<long double> > eqns, int 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 },
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why use long double when double works fine?

{ 1, 2, 3, 4 },
{ 3, -4, 0, 10 }};
EquationSolver solver;
vector<long double> ans;
try
{
ans = solver.solveByGaussJordan(equations, varc);
}
catch(runtime_error &e)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not a big fan of try catch.

{
cout << "Error found: " << e.what() << endl;
return -1;
}

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Print the matrix after the eliminations.

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 (by Backsubstitution)," << endl
<< "x == " << ans[0] << endl
<< "y == " << ans[1] << endl
<< "z == " << ans[2] << endl;
}
19 changes: 16 additions & 3 deletions contents/gaussian_elimination/gaussian_elimination.md
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -314,6 +314,9 @@ In code, this process might look like this:
{% 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)
{% sample lang="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)
Expand Down Expand Up @@ -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}
Expand All @@ -384,6 +387,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_cpp"](code/c/gaussian_elimination.c)
{% sample lang="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)
Expand All @@ -403,6 +408,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_cpp"](code/c/gaussian_elimination.c)
{% sample lang="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" %}
Expand Down Expand Up @@ -440,6 +447,8 @@ Here it is in code:
[import:67-93, lang:"julia"](code/julia/gaussian_elimination.jl)
{% sample lang="c" %}
[import:64-82, lang:"c_cpp"](code/c/gaussian_elimination.c)
{% sample lang="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)
Expand Down Expand Up @@ -481,6 +490,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_cpp"](code/c/gaussian_elimination.c)
{% sample lang="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" %}
Expand All @@ -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$$.

Expand Down Expand Up @@ -536,6 +547,8 @@ There are also plenty of other solvers that do similar things that we will get t
[import, lang:"julia"](code/julia/gaussian_elimination.jl)
{% sample lang="c" %}
[import, lang:"c_cpp"](code/c/gaussian_elimination.c)
{% sample lang="cpp" %}
[import, lang:"c_cpp"](code/c++/gaussian_elimination.cpp)
{% sample lang="rs" %}
[import, lang:"rust"](code/rust/gaussian_elimination.rs)
{% sample lang="hs" %}
Expand Down