/*
 * eliminacja gaussa v1, basic mode (musi istnieć jedno rozwiązanie)
 * specjalnie dla drg
 * (c) 2006 by tear
 * MDJCK
 */
 

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

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


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

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

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");
	}	
}

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

	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;
/*
  2    1    3    5 
  4    3    1    6 
  3    4    5    2 
*/
	dump();

	for (i = 0; i < rows; i++) {
		double t1;
		
		t1 = m(i, i);

		if (!t1) {
			abort(); /* case kiedy bedziemy mieli rozwiazanie zalezne od ilustam parametrow */
		}
		
		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;
}