/*
 * eliminacja gaussa v2, nadal nie obsluguje tylko jednorozwiazaniowe uklady
 * specjalnie dla drg
 * (c) 2006 by tear
 * MDJCK
 */
 

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

static double * matrix;
static int cols;
static int rows;


static double m(int r, int c)
{
	return matrix[r * cols + c];
}

static double * pm(int r, int c)
{
	return matrix + r * cols + c;
}

static void dump()
{
	int i, j;
	
	printf("---\n");
	for (i = 0; i < rows; i++) {
		for (j = 0; j < cols; j++) {
			printf("%.2f ", m(i, j));
		}
		printf("\n");
	}	
}

static int find_row(int n)
{
	int i;

	for (i = n; i < rows; i++) {
		if (m(i, n))
			break;
	}
	if (i >= rows)
		i = -1;
	if (i != n) {
		printf("ha!\n");
	}
	return i;
}


static void swap_rows(int a, int b)
{
	double x;
	int i;
	
	/* disclaimer: verylame */
	for (i = 0; i < cols; i++) {
		x = m(a, i);
		*pm(a, i) = m(b, i);
		*pm(b, i) = x;
	}
}

int main()
{
	int i, j, k;
	double * xp;

#ifdef DATASET1
	if (matrix)
		abort();
	cols = 4;
	rows = 3;
	matrix = malloc(cols * rows * sizeof(matrix[0]));
	xp = matrix;
	*xp++ = 2;
	*xp++ = 1;
	*xp++ = 3;
	*xp++ = 5;
	*xp++ = 4;
	*xp++ = 3;
	*xp++ = 1;
	*xp++ = 6;
	*xp++ = 3;
	*xp++ = 4;
	*xp++ = 5;
	*xp = 2;
#endif
#ifdef DATASET2
	if (matrix)
		abort();
	cols = 4;
	rows = 3;
	matrix = malloc(cols * rows * sizeof(matrix[0]));
	xp = matrix;
	*xp++ =  0; *xp++ =  6; *xp++ =  2; *xp++ = 4;
	*xp++ = -3; *xp++ = -6; *xp++ = -3; *xp++ = 5;
	*xp++ =  6; *xp++ =  3; *xp++ =  4; *xp++ = -9;
#endif
	if (!matrix)
		abort();
/*
  2    1    3    5 
  4    3    1    6 
  3    4    5    2 
*/
	dump();

	for (i = 0; i < rows; i++) {
		double t1;
		int t4;

		t4 = find_row(i);
		if (t4 == -1) {
			abort();
		}

		if (t4 != i) {
			swap_rows(t4, i);
		}

		t1 = m(i, i);

		for (j = i; j < cols; j++) {
			*pm(i, j) /= t1;
		}
		dump();
#if 1
		for (j = i + 1; j < rows; j++) {
			double t2;

			t2 = m(j, i);
			for (k = i; k < cols; k++) {
				*pm(j, k) -= t2 * m(i, k);
			}
		}
#endif
	}
	printf("---\n");
	printf("mid\n");
	for (i = rows - 1; i > 0; i--) {
		for (j = 0; j < i; j++) {
			double t3;
			
			t3 = m(j, i);

#if 0			
			if (!t3) {
				abort();
			}
#endif
			for (k = i; k < cols; k++) {
				*pm(j, k) -= t3 * m(i, k);
			}
		}
		dump();
	}
	
	
	return 0;
}