Created: 2019-12-11 ᐱ 15:32
\[ \begin{array}{l} a_1x + b_1y + c_1z = d_1 \\ a_2x + b_2y + c_2z = d_2 \\ a_3x + b_3y + c_3z = d_3 \\ \end{array} \]
1: #define MAX_N 100
2: struct AugmentedMatrix { double mat[MAX_N][MAX_N + 1]; };
3: struct ColumnVector { double vec[MAX_N]; };
4: ColumnVector GaussianElimination(int N, AugmentedMatrix Aug) {// O(N^3)
5: // input: N, Augmented Matrix Aug, output: Column vector X, the answer
6: int i, j, k, l; double t; ColumnVector X;
7: // the forward elimination phase
8: for (j = 0; j < N - 1; j++) {
9: l = j;
10: for (i = j + 1; i < N; i++) // which row has largest column value
11: if (fabs(Aug.mat[i][j]) > fabs(Aug.mat[l][j])) l = i;
12: // swap this pivot row, reason: to minimize floating point error
13: for (k = j; k <= N; k++)
14: t = Aug.mat[j][k], Aug.mat[j][k] = Aug.mat[l][k], Aug.mat[l][k] = t;
15: // the actual forward elimination phase
16: for (i = j + 1; i < N; i++)
17: for (k = N; k >= j; k--)
18: Aug.mat[i][k] -= Aug.mat[j][k] * Aug.mat[i][j] / Aug.mat[j][j];
19: }
20: for (j = N - 1; j >= 0; j--) {// the back substitution phase
21: for (t = 0.0, k = j + 1; k < N; k++) t += Aug.mat[j][k] * X.vec[k];
22: X.vec[j] = (Aug.mat[j][N] - t) / Aug.mat[j][j];
23: }
24: return X;
25: }