|
| 1 | +/* Exercise 24.4 */ |
| 2 | + |
| 3 | +#include<iostream> |
| 4 | +#include<sstream> |
| 5 | +#include<iomanip> |
| 6 | +#include"../Matrix/Matrix11.h" |
| 7 | +#include"../Matrix/MatrixIO11.h" |
| 8 | +#include"C24_Exercise_24.4.h" |
| 9 | + |
| 10 | +using namespace std; |
| 11 | + |
| 12 | +inline void error(string s) { throw runtime_error(s); } |
| 13 | +inline void error(const string& s, const string& s2) { error(s + s2); } |
| 14 | +inline void error(const string& s, int i) { ostringstream os; os << s << ": " << i; error(os.str()); } |
| 15 | +inline void keep_window_open() { char ch; cin >> ch; } |
| 16 | + |
| 17 | +int main() |
| 18 | +{ |
| 19 | + enum Action { |
| 20 | + EXIT = -1, PRINTACTIONLIST, |
| 21 | + CASE1, CASE2, CASE3, CASE4, CASE5 |
| 22 | + }; |
| 23 | + const string actionList = "\tList of actions:\n" |
| 24 | + " (1) Gaussian ilimination\n" |
| 25 | + " (-1) Exit, (0) Print the list of actions\n"; |
| 26 | + cout << actionList; |
| 27 | + int action; |
| 28 | + bool cond{ true }; |
| 29 | + while (cond) try { |
| 30 | + cout << "\nPlease enter the action: "; |
| 31 | + if (!(cin >> action)) { ClearInput(cin); error("Error. Incorrect input"); } |
| 32 | + char ch; |
| 33 | + cin.get(ch); |
| 34 | + switch (action) { |
| 35 | + case CASE1: { |
| 36 | + try { |
| 37 | + cout << endl; |
| 38 | + double a[3][3] = |
| 39 | + { |
| 40 | + {1,2,3}, |
| 41 | + {2,3,4}, |
| 42 | + {3,4,1} |
| 43 | + }; |
| 44 | + double b[3] = { 14,20,14 }; |
| 45 | + |
| 46 | + Matrix A(a); |
| 47 | + Vector B(b); |
| 48 | + |
| 49 | + cout << "Solving A * x = B" << endl; |
| 50 | + cout << "A =\n" << A << endl; |
| 51 | + cout << "B = " << B << endl; |
| 52 | + |
| 53 | + Vector x = classical_gaussian_elimination(A, B); |
| 54 | + cout << "x = " << x << endl; |
| 55 | + cout << vsp_2; |
| 56 | + |
| 57 | + |
| 58 | + double c[3][3] = |
| 59 | + { |
| 60 | + {2,1,-1}, |
| 61 | + {-3,-1,2}, |
| 62 | + {-2,1,2} |
| 63 | + }; |
| 64 | + double d[3] = { 8,-11,-3 }; |
| 65 | + |
| 66 | + Matrix C(c); |
| 67 | + Vector D(d); |
| 68 | + |
| 69 | + cout << "Solving C * y = D" << endl; |
| 70 | + cout << "C =\n" << C << endl; |
| 71 | + cout << "D = " << D << endl; |
| 72 | + |
| 73 | + Vector y = classical_gaussian_elimination(C, D); |
| 74 | + cout << "y = " << y << endl; |
| 75 | + cout << vsp_2; |
| 76 | + |
| 77 | + |
| 78 | + double e[3][3] = |
| 79 | + { |
| 80 | + {1,3,5}, |
| 81 | + {1,-2,4}, |
| 82 | + {0,1,3} |
| 83 | + }; |
| 84 | + double f[3] = { 6,-8,0 }; |
| 85 | + |
| 86 | + Matrix E(e); |
| 87 | + Vector F(f); |
| 88 | + |
| 89 | + cout << "Solving E * z = F" << endl; |
| 90 | + cout << "E =\n" << E << endl; |
| 91 | + cout << "F = " << F << endl; |
| 92 | + |
| 93 | + Vector z = classical_gaussian_elimination(E, F); |
| 94 | + cout << "z = " << z << endl; |
| 95 | + |
| 96 | + |
| 97 | + cout << vsp_2; |
| 98 | + } |
| 99 | + catch (Elim_failure& e) { |
| 100 | + cerr << e.what() << endl; |
| 101 | + } |
| 102 | + catch (Back_subst_failure& e) { |
| 103 | + cerr << e.what() << endl; |
| 104 | + } |
| 105 | + break; |
| 106 | + } |
| 107 | + case PRINTACTIONLIST: |
| 108 | + cout << actionList; |
| 109 | + break; |
| 110 | + case EXIT: |
| 111 | + cond = false; |
| 112 | + break; |
| 113 | + default: |
| 114 | + error("Error. Incorrect action number"); |
| 115 | + break; |
| 116 | + } |
| 117 | + } |
| 118 | + catch (runtime_error& e) { |
| 119 | + cerr << e.what() << endl; |
| 120 | + } |
| 121 | + catch (...) { |
| 122 | + cerr << "Error. Exception\n"; |
| 123 | + return 1; |
| 124 | + } |
| 125 | + return 0; |
| 126 | +} |
| 127 | + |
| 128 | +void ClearInput(istream& is) |
| 129 | +{ |
| 130 | + is.clear(); |
| 131 | + is.ignore(numeric_limits<streamsize>::max(), '\n'); |
| 132 | +} |
| 133 | + |
| 134 | +Vector classical_gaussian_elimination(Matrix A, Vector b) |
| 135 | +{ |
| 136 | + classical_elimination(A, b); |
| 137 | + return back_substitution(A, b); |
| 138 | +} |
| 139 | + |
| 140 | +void classical_elimination(Matrix& A, Vector& b) |
| 141 | +{ |
| 142 | + const Numeric_lib::Index n = A.dim1(); |
| 143 | + // traverse from 1st column to the next-to-last |
| 144 | + // filling zeros into all elements under the diagonal: |
| 145 | + for (Numeric_lib::Index j = 0; j < n - 1; ++j) { |
| 146 | + const double pivot = A(j, j); |
| 147 | + if (pivot == 0) throw Elim_failure("Elimination failure in row " + to_string(j)); |
| 148 | + // fill zeros into each element under the diagonal of the ith row: |
| 149 | + for (Numeric_lib::Index i = j + 1; i < n; ++i) { |
| 150 | + const double mult = A(i, j) / pivot; |
| 151 | + A[i].slice(j) = Numeric_lib::scale_and_add(A[j].slice(j), -mult, A[i].slice(j)); |
| 152 | + b(i) -= mult * b(j); // make the corresponding change to b |
| 153 | + } |
| 154 | + } |
| 155 | +} |
| 156 | + |
| 157 | +Vector back_substitution(const Matrix& A, const Vector& b) |
| 158 | +{ |
| 159 | + const Numeric_lib::Index n = A.dim1(); |
| 160 | + Vector x(n); |
| 161 | + for (Numeric_lib::Index i = n - 1; i >= 0; --i) { |
| 162 | + double s = b(i) - Numeric_lib::dot_product(A[i].slice(i + 1), x.slice(i + 1)); |
| 163 | + if (double m = A(i, i)) x(i) = s / m; |
| 164 | + else throw Back_subst_failure("Back substitution failure in row " + to_string(i)); |
| 165 | + } |
| 166 | + return x; |
| 167 | +} |
0 commit comments