Last active
October 3, 2025 18:56
-
-
Save kevyonan/c079e446d96e6d7d830061d202ff0259 to your computer and use it in GitHub Desktop.
testing gaussian elimination to solve system of equations.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #include <stdio.h> | |
| #include <math.h> | |
| size_t idx_2_to_1(size_t const idx_dim1, size_t const idx_dim2, size_t const dim_size) { | |
| return (idx_dim1 * dim_size) + idx_dim2; | |
| } | |
| enum RREFResult { | |
| RREFResultOk, | |
| RREFResultFreeVars, | |
| RREFResultBadMatrix, | |
| }; | |
| enum RREFResult gaussian_rref( | |
| size_t const matrix_len, | |
| double A[const restrict static matrix_len*matrix_len], | |
| double v[const restrict static matrix_len] | |
| ) { | |
| size_t curr_pivot_row = 0; | |
| for( size_t c=0; c < matrix_len && curr_pivot_row < matrix_len; c++ ) { | |
| size_t curr_pivot = curr_pivot_row; | |
| size_t const rc = idx_2_to_1(curr_pivot_row, c, matrix_len); | |
| double cur_max = fabs(A[rc]); | |
| for( size_t i = curr_pivot_row+1; i < matrix_len; i++ ) { | |
| size_t const ic = idx_2_to_1(i, c, matrix_len); | |
| double const mag = fabs(A[ic]); | |
| if( cur_max < mag ) { | |
| cur_max = mag; curr_pivot = i; | |
| } | |
| } | |
| if( cur_max<1e-15 ) { | |
| continue; | |
| } | |
| if( curr_pivot != curr_pivot_row ) { | |
| double const tmpv = v[curr_pivot_row]; | |
| v[curr_pivot_row] = v[curr_pivot]; | |
| v[curr_pivot] = tmpv; | |
| for( size_t j=0; j < matrix_len; j++ ) { | |
| size_t const rj = idx_2_to_1(curr_pivot_row, j, matrix_len); | |
| size_t const pj = idx_2_to_1(curr_pivot, j, matrix_len); | |
| double const tmpA = A[rj]; | |
| A[rj] = A[pj]; | |
| A[pj] = tmpA; | |
| } | |
| } | |
| size_t const rc2 = idx_2_to_1(curr_pivot_row, c, matrix_len); | |
| double const piv = A[rc2]; | |
| for( size_t j=0; j < matrix_len; j++ ) { | |
| size_t const rj = idx_2_to_1(curr_pivot_row, j, matrix_len); | |
| A[rj] = (A[rj] / piv); | |
| } | |
| v[curr_pivot_row] = (v[curr_pivot_row] / piv); | |
| for( size_t i=0; i < matrix_len; i++ ) { | |
| if( i==curr_pivot_row ) { | |
| continue; | |
| } | |
| size_t const ic = idx_2_to_1(i, c, matrix_len); | |
| double const factor = A[ic]; | |
| /// If already ~0, skip the row update. | |
| if( fabs(factor) < 1e-15 ) { | |
| continue; | |
| } | |
| for( size_t j=0; j < matrix_len; j++ ) { | |
| size_t const ij = idx_2_to_1(i, j, matrix_len); | |
| size_t const rj = idx_2_to_1(curr_pivot_row, j, matrix_len); | |
| A[ij] = (A[ij] - (factor * A[rj])); | |
| } | |
| v[i] = (v[i] - (factor * v[curr_pivot_row])); | |
| /// A[i,c] is now (numerically) zero. | |
| } | |
| curr_pivot_row++; | |
| } | |
| // 6) Detect inconsistency: a zero row of A with nonzero RHS. | |
| for( size_t i=0; i < matrix_len; i++ ) { | |
| bool all_zero = true; | |
| for( size_t j=0; j < matrix_len; j++ ) { | |
| double const aij = A[idx_2_to_1(i, j, matrix_len)]; | |
| if( !(fabs(aij) < 1e-15) ) { | |
| all_zero = false; | |
| break; | |
| } | |
| } | |
| if( all_zero && !(fabs(v[i]) < 1e-15) ) { | |
| return RREFResultBadMatrix; /// inconsistent | |
| } | |
| } | |
| return curr_pivot_row==matrix_len? RREFResultOk : RREFResultFreeVars; | |
| } | |
| int main() { | |
| enum { | |
| M_SIZE = 3, | |
| }; | |
| /// [1 4 6, 2 3 4, 5 2 1] | |
| double matrix[M_SIZE][M_SIZE] = { | |
| {-6,1,1}, | |
| {1,-2,2}, | |
| {1,1,-1}, | |
| }; | |
| double v[M_SIZE] = {11, 4, 1}; | |
| printf("Matrix :: \n"); | |
| for( size_t row=0; row < M_SIZE; row++ ) { | |
| for( size_t col=0; col < M_SIZE; col++ ) { | |
| printf("| %+f |", matrix[row][col]); | |
| } | |
| printf("| %+f\n", v[row]); | |
| } | |
| gaussian(M_SIZE, matrix, v); | |
| printf("\nSolution :: \n[x]: %+f\n[y]: %+f\n[z]: %+f\n", v[0], v[1], v[2]); | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment