gistfile1.c

/*
 ============================================================================
 Name        : gausstest.c
 Author      : 
 Version     :
 Copyright   : Your copyright notice
 Description : Hello World in C, Ansi-style
 ============================================================================
 */

#include <stdio.h>
#include <stdlib.h>

//  swaprows - exchanges the contents of row0 and row1 in a 2d array
void swaprows(double** arr, long row0, long row1) {
	double* temp;
	temp = arr[row0];
	arr[row0] = arr[row1];
	arr[row1] = temp;
}

/**
 * WIP
 *
 * Standard gaussian elimination with partial pivoting
 *
 * NB: Destructive!
 */
void gelim(double** arr, long nrows, long ncols) {

	//number of forward eliminations = ncols - 1
	long col;
	for (col = 0; col < ncols-1; ++col) {

		long pivotrow = col;

		//if matrix too short, we must abort
		if (pivotrow >= nrows) {
			return;
		}

		// following is the pivoting.
		// pivoting is done to improve the numerical stability of the algorithm,
		// which is a fancy way of saying we want to reduce rounding errors.
		// the basic idea here is that as eliminator we use the biggest coefficient
		// in the column. The following search-and-swap is essentially equivalent to
		// sorting rows by coefficients, descending, sorting by the first column first,
		// then the second, and so on.

		// initialize maxer
		double max = arr[pivotrow][col];
		long maxrow = pivotrow;
		long row;
		// find biggest value in remaining matrix
		for (row = pivotrow+1; row < nrows; ++row) {
			if (abs(arr[row][col]) > abs(max)) {
				max = arr[row][col];
				maxrow = row;
			}
		}

		//swap max row up to pivot row
		if (maxrow != pivotrow) {
			swaprows(arr,maxrow,pivotrow);
		}

		//eliminate in all rows under the pivot row
		for(row=pivotrow+1;row<nrows;++row) {

			// following is the crucial elimination step.
			// what we do here is that we want to zero out the
			// coefficient at arr[row][col]. Because the only action we are allowed
			// is to add or subtract a multiple of any row to/from another row,
			// we have to take our pivot row and multiply it by a factor that
			// will make it so the coefficient in the eliminating row equals the
			// coefficient in the to-be-eliminated row.
			// the solution of arr[pivotrow][col] * x = arr[row][col] is of course
			// x = arr[row][col] / arr[pivotrow][col], and that is our factor.

			// if arr[pivotrow][col] is zero, all coefficients in this column are zero.

			if (0 == arr[pivotrow][col]) {
				// http://dl.dropbox.com/u/2808338/nothing-to-do-here-template.jpg.scaled500.jpg
				continue;
			}

			double factor = arr[row][col] / arr[pivotrow][col];


			long elimcol;
			for (elimcol=0;elimcol<ncols;++elimcol) {
				arr[row][elimcol] = arr[row][elimcol] -  arr[pivotrow][elimcol] * factor;
			}
		}
	}
}

int main(void) {

	//allocation
	long nrows = 3;
	long ncols = 3;

	double** arr2 = (double**) malloc(nrows * sizeof(double*));
	long row;
	for (row = 0; row < nrows; ++row) {
		arr2[row] = (double*) malloc(ncols * sizeof(double));
	}

	arr2[0][0] = 3; arr2[0][1] = 3; arr2[0][2] = 3;
	arr2[1][0] = 3; arr2[1][1] = 3; arr2[1][2] = 3;
	arr2[2][0] = 3; arr2[2][1] = 3; arr2[2][2] = 3;

	printf("%8.3f %8.3f %8.3f\n",arr2[0][0],arr2[0][1],arr2[0][2]);
	printf("%8.3f %8.3f %8.3f\n",arr2[1][0],arr2[1][1],arr2[1][2]);
	printf("%8.3f %8.3f %8.3f\n",arr2[2][0],arr2[2][1],arr2[2][2]);

	gelim(arr2,nrows,ncols);

	printf("\n");

	printf("%8.3f %8.3f %8.3f\n",arr2[0][0],arr2[0][1],arr2[0][2]);
	printf("%8.3f %8.3f %8.3f\n",arr2[1][0],arr2[1][1],arr2[1][2]);
	printf("%8.3f %8.3f %8.3f\n",arr2[2][0],arr2[2][1],arr2[2][2]);


	return EXIT_SUCCESS;
}