Study of the algebra of vector spaces.

Gaussian Elimination

It is a fairly simple thing to implement the algorithm for Gaussian Elimination in a manner that produces reasonably stable results. Perfect for reducing simple matrices.

#include <stdio.h>
#include <stdint.h>
#include <math.h>

typedef double f64;

void MatrixPrint(f64* matrix, int row_count, int col_count) {
	for (int row_idx = 0; row_idx < row_count; ++row_idx) {
		printf("[");
		for (int col_idx = 0; col_idx < col_count; ++col_idx) {
			if (col_idx > 0) {
				printf("\t");
			}
			printf("%0.02f", matrix[col_idx + (row_idx * col_count)]);
		}
		printf("]\n");
	}
}

void MatrixSwapRows(f64* matrix, int row_count, int col_count, int swap_row_a, int swap_row_b) {
	f64* row_a = &matrix[swap_row_a * col_count];
	f64* row_b = &matrix[swap_row_b * col_count];

	for (int col_idx = 0; col_idx < col_count; ++col_idx) {
		f64 tmp = row_a[col_idx];
		row_a[col_idx] = row_b[col_idx];
		row_b[col_idx] = tmp;
	}
}

void MatrixRowReduce(f64* matrix, int row_count, int col_count) {
#define MatE(R,C) matrix[(C) + ((R) * col_count)]
	/* Reduce the matrix to row-echelon form */
	int pivot_row = 0;
	int pivot_col = 0;
  
	while (pivot_row < row_count && pivot_col < col_count) {
		/* Find the row with the largest value in its pivot_col to use as next pivot_row */
		int next_row_idx = -1;
		{
			f64 last_value = 0;
			for (int row_idx = pivot_row; row_idx < row_count; ++row_idx) {
				f64 value = MatE(row_idx, pivot_col);
        
				if (fabs(value) == 0.0) {
					continue;
				}
        
				if (fabs(value) > fabs(last_value)) {
					last_value = value;
					next_row_idx = row_idx;
				}
			}
		}
    
		if (next_row_idx == -1) {
			/* No pivot in this column so pass to the next column */
			++pivot_col;
		} else {
			MatrixSwapRows(matrix, row_count, col_count, pivot_row, next_row_idx);
      
			/* Eliminate all rows below the pivot */
			for (int row_idx = pivot_row + 1; row_idx < row_count; ++row_idx) {
				f64 ratio = MatE(row_idx, pivot_col) / MatE(pivot_row, pivot_col);
				/* Clear current element and update all the elements after it in this row */
				MatE(row_idx, pivot_col) = 0.0;
				for (int col_idx = pivot_col + 1; col_idx < col_count; ++col_idx) {
					MatE(row_idx, col_idx) -= ratio * MatE(pivot_row, col_idx);
				}
			}
      
			++pivot_row;
			++pivot_col;
		}
	}
  
	/* Eliminate upwards */
	for (int row_idx = 1; row_idx < row_count; ++row_idx) {
		int curr_col_idx = row_idx;
		f64 ratio = MatE(row_idx - 1, curr_col_idx) / MatE(row_idx, curr_col_idx);
		matrix[curr_col_idx + ((row_idx - 1) * col_count)] = 0.0;
		for (int col_idx = curr_col_idx + 1; col_idx < col_count; ++col_idx) {
			MatE(row_idx - 1, col_idx) -= ratio * MatE(row_idx, col_idx);
		}
	}
  
	/* Finally normalize the rows so we have it in row-reduced echelon form */
	for (int row_idx = 0; row_idx < row_count; ++row_idx) {
		f64 value = matrix[row_idx + (row_idx * col_count)];
		if (fabs(value) == 0.0) {
			continue;
		}
    
		f64 ratio = 1.0 / value;
		for (int col_idx = row_idx; col_idx < col_count; ++col_idx) {
			matrix[col_idx + (row_idx * col_count)] *= ratio;
		}
	}

#undef MatE
}

incoming(1): oscillators