From 2f52647a0e172e34fd5351c4d2912180697cb735 Mon Sep 17 00:00:00 2001 From: Anton Te Date: Mon, 16 Jul 2018 00:26:59 -0700 Subject: [PATCH 1/3] Implement Gaussian Elimination in C++ --- .../code/c++/gaussian_elimination.cpp | 133 ++++++++++++++++++ .../gaussian_elimination.md | 8 ++ 2 files changed, 141 insertions(+) create mode 100644 contents/gaussian_elimination/code/c++/gaussian_elimination.cpp diff --git a/contents/gaussian_elimination/code/c++/gaussian_elimination.cpp b/contents/gaussian_elimination/code/c++/gaussian_elimination.cpp new file mode 100644 index 000000000..fff28e2b3 --- /dev/null +++ b/contents/gaussian_elimination/code/c++/gaussian_elimination.cpp @@ -0,0 +1,133 @@ +#include +#include +#include +#include +#include + +void gaussian_elimination(std::vector& a, int cols) { + assert(a.size() % cols == 0); + int rows = a.size() / cols; + + int row = 0; + + // Main loop going through all columns + for (int col = 0; col < cols - 1; ++col) { + // Step 1: finding the maximum element for each column + int max_index = [&]() { + int res = row; + int max_element = a[col + row * cols]; + for (int r = row + 1; r < rows; ++r) + if (max_element < std::abs(a[col + r * cols])) { + max_element = std::abs(a[col + r * cols]); + res = r; + } + return res; + }(); + + // Check to make sure matrix is good! + if (a[col + max_index * cols] == 0) { + std::cout << "matrix is singular!\n"; + continue; + } + + // Step 2: swap row with highest value for that column to the top + for (int c = 0; c < cols; ++c) + std::swap(a[c + row * cols], a[c + max_index * cols]); + + // Loop for all remaining rows + for (int i = row + 1; i < rows; ++i) { + + // Step 3: finding fraction + auto fraction = a[col + i * cols] / a[col + row * cols]; + + // loop through all columns for that row + for (int j = col + 1; j < cols; ++j) { + + // Step 4: re-evaluate each element + a[j + i * cols] -= a[j + row * cols] * fraction; + } + + // Step 5: Set lower elements to 0 + a[col + i * cols] = 0; + } + ++row; + } +} + +std::vector back_substitution(const std::vector& a, int cols) { + assert(a.size() % cols == 0); + int rows = a.size() / cols; + + // Creating the solution Vector + std::vector soln(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; +} + +void gauss_jordan_elimination(std::vector& a, int cols) { + assert(a.size() % cols == 0); + // After this, we know what row to start on (r-1) + // to go back through the matrix + int row = 0; + for (int col = 0; col < cols - 1; ++col) { + if (a[col + row * cols] != 0) { + + // divide row by pivot and leaving pivot as 1 + for (int i = cols - 1; i >= static_cast(col); --i) + a[i + row * cols] /= a[col + row * cols]; + + // subtract value from above row and set values above pivot to 0 + for (int i = 0; i < static_cast(row - 1); ++i) + for (int j = cols - 1; j >= static_cast(col); --j) + a[j + i * cols] -= a[col + i * cols] * a[j + row * cols]; + ++row; + } + } +} + +void print_matrix(const std::vector& a, int cols) { + assert(a.size() % cols == 0); + int rows = a.size() / cols; + for (int i = 0; i < rows; ++i) { + std::cout << "["; + for (int j = 0; j < cols; ++j) { + std::cout << a[j + i * cols] << " "; + } + std::cout << "]\n"; + } +} + +int main() { + std::vector a = {2, 3, 4, 6, 1, 2, 3, 4, 3, -4, 0, 10}; + const int cols = 4; + assert(a.size() % cols == 0); + + gaussian_elimination(a, cols); + print_matrix(a, cols); + + auto soln = back_substitution(a, 4); + + for (auto element : soln) + std::cout << element << std::endl; + + gauss_jordan_elimination(a, cols); + print_matrix(a, cols); + + soln = back_substitution(a, 4); + + for (auto element : soln) + std::cout << element << std::endl; +} diff --git a/contents/gaussian_elimination/gaussian_elimination.md b/contents/gaussian_elimination/gaussian_elimination.md index d3c3853c9..ab36a2ce9 100644 --- a/contents/gaussian_elimination/gaussian_elimination.md +++ b/contents/gaussian_elimination/gaussian_elimination.md @@ -356,6 +356,8 @@ In code, this looks like: [import:1-45, 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-54, lang:"c_cpp"](code/c++/gaussian_elimination.cpp) {% sample lang="rs" %} [import:41-78, lang:"rust"](code/rust/gaussian_elimination.rs) {% sample lang="hs" %} @@ -389,6 +391,8 @@ Here it is in code: {% sample lang="c" %} This code does not exist yet in C, so here's Julia code (sorry for the inconvenience) [import:67-93, lang:"julia"](code/julia/gaussian_elimination.jl) +{% sample lang="cpp" %} +[import:79-98, 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) @@ -420,6 +424,8 @@ In code, this involves keeping a rolling sum of all the values we substitute in [import:47-67, lang:"julia"](code/julia/gaussian_elimination.jl) {% sample lang="c" %} [import:46-58, lang:"c_cpp"](code/c/gaussian_elimination.c) +{% sample lang="cpp" %} +[import:56-77, lang:"c_cpp"](code/c++/gaussian_elimination.cpp) {% sample lang="rs" %} [import:79-94, lang:"rust"](code/rust/gaussian_elimination.rs) {% sample lang="hs" %} @@ -442,6 +448,8 @@ As for what's next... Well, we are in for a treat! The above algorithm clearly h [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" %} From e84341861eaa58618a0388c76677ef4033ee00a1 Mon Sep 17 00:00:00 2001 From: Anton Te Date: Sat, 21 Jul 2018 14:04:13 -0700 Subject: [PATCH 2/3] Address code review comments: - add curly bracket - remove static_cast - remove IIFE - fix format for the matrix - format output with tabs - properly sanitize user input - remove extraneous assignment - fix off-by-1 --- .../code/c++/gaussian_elimination.cpp | 39 +++++++++---------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/contents/gaussian_elimination/code/c++/gaussian_elimination.cpp b/contents/gaussian_elimination/code/c++/gaussian_elimination.cpp index fff28e2b3..745e13356 100644 --- a/contents/gaussian_elimination/code/c++/gaussian_elimination.cpp +++ b/contents/gaussian_elimination/code/c++/gaussian_elimination.cpp @@ -3,6 +3,7 @@ #include #include #include +#include void gaussian_elimination(std::vector& a, int cols) { assert(a.size() % cols == 0); @@ -13,16 +14,12 @@ void gaussian_elimination(std::vector& a, int cols) { // Main loop going through all columns for (int col = 0; col < cols - 1; ++col) { // Step 1: finding the maximum element for each column - int max_index = [&]() { - int res = row; - int max_element = a[col + row * cols]; - for (int r = row + 1; r < rows; ++r) - if (max_element < std::abs(a[col + r * cols])) { - max_element = std::abs(a[col + r * cols]); - res = r; - } - return res; - }(); + int max_index = row; + for (int r = row + 1; r < rows; ++r) { + if (a[col + max_index * cols] < std::abs(a[col + r * cols])) { + max_index = r; + } + } // Check to make sure matrix is good! if (a[col + max_index * cols] == 0) { @@ -61,10 +58,6 @@ std::vector back_substitution(const std::vector& a, int cols) { // Creating the solution Vector std::vector soln(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) { @@ -86,12 +79,12 @@ void gauss_jordan_elimination(std::vector& a, int cols) { if (a[col + row * cols] != 0) { // divide row by pivot and leaving pivot as 1 - for (int i = cols - 1; i >= static_cast(col); --i) + for (int i = cols - 1; i >= col; --i) a[i + row * cols] /= a[col + row * cols]; // subtract value from above row and set values above pivot to 0 - for (int i = 0; i < static_cast(row - 1); ++i) - for (int j = cols - 1; j >= static_cast(col); --j) + for (int i = 0; i < row; ++i) + for (int j = cols - 1; j >= col; --j) a[j + i * cols] -= a[col + i * cols] * a[j + row * cols]; ++row; } @@ -104,16 +97,22 @@ void print_matrix(const std::vector& a, int cols) { for (int i = 0; i < rows; ++i) { std::cout << "["; for (int j = 0; j < cols; ++j) { - std::cout << a[j + i * cols] << " "; + std::cout << std::fixed << a[j + i * cols] << "\t"; } std::cout << "]\n"; } } int main() { - std::vector a = {2, 3, 4, 6, 1, 2, 3, 4, 3, -4, 0, 10}; + std::vector a = { 2, 3, 4, 6, + 1, 2, 3, 4, + 3, -4, 0, 10 }; const int cols = 4; - assert(a.size() % cols == 0); + if (a.size() % cols != 0) + { + std::cout << "The input dimentions are incorrect\n"; + return 1; + } gaussian_elimination(a, cols); print_matrix(a, cols); From dc87fe0d2c10866c7080f5c193e10d52a8bf9b15 Mon Sep 17 00:00:00 2001 From: Anton Te Date: Sat, 4 Aug 2018 21:59:44 -0700 Subject: [PATCH 3/3] Address code review comments - remove asserts - don't swap if it is a same row --- .../code/c++/gaussian_elimination.cpp | 10 ++++------ contents/gaussian_elimination/gaussian_elimination.md | 6 +++--- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/contents/gaussian_elimination/code/c++/gaussian_elimination.cpp b/contents/gaussian_elimination/code/c++/gaussian_elimination.cpp index 745e13356..50cc3b5a9 100644 --- a/contents/gaussian_elimination/code/c++/gaussian_elimination.cpp +++ b/contents/gaussian_elimination/code/c++/gaussian_elimination.cpp @@ -6,7 +6,6 @@ #include void gaussian_elimination(std::vector& a, int cols) { - assert(a.size() % cols == 0); int rows = a.size() / cols; int row = 0; @@ -28,8 +27,10 @@ void gaussian_elimination(std::vector& a, int cols) { } // Step 2: swap row with highest value for that column to the top - for (int c = 0; c < cols; ++c) - std::swap(a[c + row * cols], a[c + max_index * cols]); + if (row != max_index) { + for (int c = 0; c < cols; ++c) + std::swap(a[c + row * cols], a[c + max_index * cols]); + } // Loop for all remaining rows for (int i = row + 1; i < rows; ++i) { @@ -52,7 +53,6 @@ void gaussian_elimination(std::vector& a, int cols) { } std::vector back_substitution(const std::vector& a, int cols) { - assert(a.size() % cols == 0); int rows = a.size() / cols; // Creating the solution Vector @@ -71,7 +71,6 @@ std::vector back_substitution(const std::vector& a, int cols) { } void gauss_jordan_elimination(std::vector& a, int cols) { - assert(a.size() % cols == 0); // After this, we know what row to start on (r-1) // to go back through the matrix int row = 0; @@ -92,7 +91,6 @@ void gauss_jordan_elimination(std::vector& a, int cols) { } void print_matrix(const std::vector& a, int cols) { - assert(a.size() % cols == 0); int rows = a.size() / cols; for (int i = 0; i < rows; ++i) { std::cout << "["; diff --git a/contents/gaussian_elimination/gaussian_elimination.md b/contents/gaussian_elimination/gaussian_elimination.md index ab36a2ce9..bbc1acd6c 100644 --- a/contents/gaussian_elimination/gaussian_elimination.md +++ b/contents/gaussian_elimination/gaussian_elimination.md @@ -357,7 +357,7 @@ In code, this looks like: {% sample lang="c" %} [import:13-44, lang:"c_cpp"](code/c/gaussian_elimination.c) {% sample lang="cpp" %} -[import:6-54, lang:"c_cpp"](code/c++/gaussian_elimination.cpp) +[import:8-53, lang:"c_cpp"](code/c++/gaussian_elimination.cpp) {% sample lang="rs" %} [import:41-78, lang:"rust"](code/rust/gaussian_elimination.rs) {% sample lang="hs" %} @@ -392,7 +392,7 @@ Here it is in code: This code does not exist yet in C, so here's Julia code (sorry for the inconvenience) [import:67-93, lang:"julia"](code/julia/gaussian_elimination.jl) {% sample lang="cpp" %} -[import:79-98, lang:"c_cpp"](code/c++/gaussian_elimination.cpp) +[import:73-91, 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) @@ -425,7 +425,7 @@ In code, this involves keeping a rolling sum of all the values we substitute in {% sample lang="c" %} [import:46-58, lang:"c_cpp"](code/c/gaussian_elimination.c) {% sample lang="cpp" %} -[import:56-77, lang:"c_cpp"](code/c++/gaussian_elimination.cpp) +[import:55-71, lang:"c_cpp"](code/c++/gaussian_elimination.cpp) {% sample lang="rs" %} [import:79-94, lang:"rust"](code/rust/gaussian_elimination.rs) {% sample lang="hs" %}