|
| 1 | +#include <algorithm> |
| 2 | +#include <cassert> |
| 3 | +#include <cmath> |
| 4 | +#include <iostream> |
| 5 | +#include <vector> |
| 6 | + |
| 7 | +void gaussian_elimination(std::vector<double>& a, int cols) { |
| 8 | + assert(a.size() % cols == 0); |
| 9 | + int rows = a.size() / cols; |
| 10 | + |
| 11 | + int row = 0; |
| 12 | + |
| 13 | + // Main loop going through all columns |
| 14 | + for (int col = 0; col < cols - 1; ++col) { |
| 15 | + // Step 1: finding the maximum element for each column |
| 16 | + int max_index = [&]() { |
| 17 | + int res = row; |
| 18 | + int max_element = a[col + row * cols]; |
| 19 | + for (int r = row + 1; r < rows; ++r) |
| 20 | + if (max_element < std::abs(a[col + r * cols])) { |
| 21 | + max_element = std::abs(a[col + r * cols]); |
| 22 | + res = r; |
| 23 | + } |
| 24 | + return res; |
| 25 | + }(); |
| 26 | + |
| 27 | + // Check to make sure matrix is good! |
| 28 | + if (a[col + max_index * cols] == 0) { |
| 29 | + std::cout << "matrix is singular!\n"; |
| 30 | + continue; |
| 31 | + } |
| 32 | + |
| 33 | + // Step 2: swap row with highest value for that column to the top |
| 34 | + for (int c = 0; c < cols; ++c) |
| 35 | + std::swap(a[c + row * cols], a[c + max_index * cols]); |
| 36 | + |
| 37 | + // Loop for all remaining rows |
| 38 | + for (int i = row + 1; i < rows; ++i) { |
| 39 | + |
| 40 | + // Step 3: finding fraction |
| 41 | + auto fraction = a[col + i * cols] / a[col + row * cols]; |
| 42 | + |
| 43 | + // loop through all columns for that row |
| 44 | + for (int j = col + 1; j < cols; ++j) { |
| 45 | + |
| 46 | + // Step 4: re-evaluate each element |
| 47 | + a[j + i * cols] -= a[j + row * cols] * fraction; |
| 48 | + } |
| 49 | + |
| 50 | + // Step 5: Set lower elements to 0 |
| 51 | + a[col + i * cols] = 0; |
| 52 | + } |
| 53 | + ++row; |
| 54 | + } |
| 55 | +} |
| 56 | + |
| 57 | +std::vector<double> back_substitution(const std::vector<double>& a, int cols) { |
| 58 | + assert(a.size() % cols == 0); |
| 59 | + int rows = a.size() / cols; |
| 60 | + |
| 61 | + // Creating the solution Vector |
| 62 | + std::vector<double> soln(rows); |
| 63 | + |
| 64 | + // initialize the final element |
| 65 | + soln[rows - 1] = |
| 66 | + a[cols - 1 + (rows - 1) * cols] / a[cols - 1 - 1 + (rows - 1) * cols]; |
| 67 | + |
| 68 | + for (int i = rows - 1; i >= 0; --i) { |
| 69 | + auto sum = 0.0; |
| 70 | + for (int j = cols - 2; j > i; --j) { |
| 71 | + sum += soln[j] * a[j + i * cols]; |
| 72 | + } |
| 73 | + |
| 74 | + soln[i] = (a[cols - 1 + i * cols] - sum) / a[i + i * cols]; |
| 75 | + } |
| 76 | + |
| 77 | + return soln; |
| 78 | +} |
| 79 | + |
| 80 | +void gauss_jordan_elimination(std::vector<double>& a, int cols) { |
| 81 | + assert(a.size() % cols == 0); |
| 82 | + // After this, we know what row to start on (r-1) |
| 83 | + // to go back through the matrix |
| 84 | + int row = 0; |
| 85 | + for (int col = 0; col < cols - 1; ++col) { |
| 86 | + if (a[col + row * cols] != 0) { |
| 87 | + |
| 88 | + // divide row by pivot and leaving pivot as 1 |
| 89 | + for (int i = cols - 1; i >= static_cast<int>(col); --i) |
| 90 | + a[i + row * cols] /= a[col + row * cols]; |
| 91 | + |
| 92 | + // subtract value from above row and set values above pivot to 0 |
| 93 | + for (int i = 0; i < static_cast<int>(row - 1); ++i) |
| 94 | + for (int j = cols - 1; j >= static_cast<int>(col); --j) |
| 95 | + a[j + i * cols] -= a[col + i * cols] * a[j + row * cols]; |
| 96 | + ++row; |
| 97 | + } |
| 98 | + } |
| 99 | +} |
| 100 | + |
| 101 | +void print_matrix(const std::vector<double>& a, int cols) { |
| 102 | + assert(a.size() % cols == 0); |
| 103 | + int rows = a.size() / cols; |
| 104 | + for (int i = 0; i < rows; ++i) { |
| 105 | + std::cout << "["; |
| 106 | + for (int j = 0; j < cols; ++j) { |
| 107 | + std::cout << a[j + i * cols] << " "; |
| 108 | + } |
| 109 | + std::cout << "]\n"; |
| 110 | + } |
| 111 | +} |
| 112 | + |
| 113 | +int main() { |
| 114 | + std::vector<double> a = {2, 3, 4, 6, 1, 2, 3, 4, 3, -4, 0, 10}; |
| 115 | + const int cols = 4; |
| 116 | + assert(a.size() % cols == 0); |
| 117 | + |
| 118 | + gaussian_elimination(a, cols); |
| 119 | + print_matrix(a, cols); |
| 120 | + |
| 121 | + auto soln = back_substitution(a, 4); |
| 122 | + |
| 123 | + for (auto element : soln) |
| 124 | + std::cout << element << std::endl; |
| 125 | + |
| 126 | + gauss_jordan_elimination(a, cols); |
| 127 | + print_matrix(a, cols); |
| 128 | + |
| 129 | + soln = back_substitution(a, 4); |
| 130 | + |
| 131 | + for (auto element : soln) |
| 132 | + std::cout << element << std::endl; |
| 133 | +} |
0 commit comments