/*	gaussalg.c -- Matthias Jauernig, 12.12.03					*/
/*	Programm berechnet die Lösungen eines LGS mit (n,n)-Koeffizientenmatrix		*/

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

/* -- printline() -- gibt eine Linie der Länge n mit dem Zeichen z aus ----------------	*/
void printline(int l, char z){
	int i;
	for(i=0; i<=l; i++)
		printf("%c",z);
	printf("\n");
}

/* -- printhead() -- gibt die Kopfzeile der Tabelle aus -------------------------------	*/
void printhead(int n){
	int i;
	char str[10];
	printf("\n");
	printline(10*(n+1), '=');
	for(i=0; i<n; i++){
		snprintf(str, 10, "x[%d]", i);
		printf("|%8s ", str);
	}
	printf("|     b   |\n");
	printline(10*(n+1), '=');
}

/* -- printgauss() -- gibt die aktuelle Koeffizientenmatrix a in Tabellenform aus -----	*/
void printgauss(float **a, int n, char z){
	int i, j;
	for(i=0; i<n; i++){
			for(j=0; j<=n; j++)
				printf("| %7.2g ", a[i][j]);
			printf("|\n");
		}
		printline(10*(n+1), z);
}

float absol(float z){
	if(z<0.0)	return -z;
	else 		return z;
}

/* == main() ==========================================================================	*/
int main(void){
	// **a: Koeffizientenmatrix; *x: Lösungsvektor
	// quot: Quotient a[i][k]/a[k][k]; sub: für Lösungsvektor zu berechnende Summe
	// i,j,k: Laufvariablen; n: Zahl für die (n,n)-Koeffizientenmatrix
	// flag: zeigt an, ob eine Zeile vor Ende des Algorithmus' komplett 0 ist
	float **a, *x, quot, sub;
	int i, j, k, n;
	char flag;

	//################################# Eingabe ######################################
	printf(	"--------------------------------------\n"
		"| Gauss-Algorithmus ohne Pivot-Suche |\n"
		"--------------------------------------\n\n");
	do{
		printf("Wie groß ist n für die (n,n)-Koeffizientenmatrix?: ");
		scanf("%d", &n);
	}while(n<=0 && printf("!! n muss größer als 0 sein!\n\n"));

	// Speicher reservieren
	a=(float **)malloc(n*sizeof(float *));
	for(i=0; i<n; i++)
		a[i]=(float *)malloc((n+1)*sizeof(float));

	printf("\nWerte der Koeffizientenmatrix eingeben:\n");
	for(i=0; i<n; i++)
		for(j=0; j<n; j++)
			do{
				printf("Wert von a[%d][%d] eingeben: ", i, j);
				scanf("%f", &a[i][j]);
			}while(	i==j && a[i][j]==0
				&& printf("!! a[i][i] darf nicht 0 sein!\n\n"));

	printf("\nWerte des Vektors b eingeben:\n");
	for(i=0; i<n; i++){
		printf("Wert von b[%d] eingeben: ", i);
		scanf("%f", &a[i][n]);
	}

	//######################## a auf Dreiecksform bringen ############################
	if(n<=6){
		printhead(n);
		printgauss(a, n, '~');
	}

	//1. Schleife: gehe alle Zeilen von 0 bis n-1 durch, höre auf, wenn eine Zeile komplett 0
	//2. Schleife: subtrahiere von jeder Zeile unter k
	//3. Schleife: subtrahiere von jeder Spalte
	for(k=0, flag=0; k<n-1 && flag==0; k++){
		for(i=k+1; i<n; i++){
			quot=a[i][k]/a[k][k];
			for(j=k; j<=n; j++){
				a[i][j]-=a[k][j]*quot;
				if(absol(a[i][j])<1e-05) a[i][j]=0.0; //Fehler durch Rundung
			}
		}
		if(n<=6) printgauss(a, n, (k==n-2) ? '=' : '-');

		//folgender Codeteil bewirkt, dass terminiert wird, wenn eine Zeile 0 ist;
		//ist in dem Fall b[i]=0, so wird flag=1 gesetzt -> unendlich viele Lsg.
		//ist in dem Fall b[i]!=0, so wird flag=2 gesetzt -> keine Lösungen
		for(i=k+1; i<n && flag!=2; i++){
			for(j=k, flag=1; j<n && flag==1; j++)
				if(a[i][j]!=0) flag=0;
			if(flag==1)
				if(a[i][n]!=0)
					flag=2;
		}
		if(flag!=0){
			if(flag==1)
				printf("=> %d. Zeile 0=0: also unendlich viele Lösungen!\n", i);
			else
				printf("=> %d. Zeile 0!=%g: also keine Lösung!\n", i, a[i-1][n]);
		}

	}

	//################### durch Rücksubstitution x[i] berechnen ######################
	if(flag==0){
		x=(float *)malloc(n*sizeof(float));
		printf("\nLösungen nach dem Gaussschen Algorithmus:\n");
		for(i=n-1; i>=0; i--){
			sub=0;
			//Summe berechnen, die abgezogen wird
			for(j=i+1; j<n; j++)
				sub+=a[i][j]*x[j];
			x[i]=(a[i][n]-sub)/a[i][i];
			printf("=> x[%d] = %g\n", i, x[i]);
		}
		free(x);
	}
	printf("\n\n");

	//################################# Aufräumen ####################################
	for(i=0; i<n; i++)
		free(a[i]);
	free(a);
	return 0;
}

/* Beispielabläufe
##########################################################################################
--------------------------------------
| Gauss-Algorithmus ohne Pivot-Suche |
--------------------------------------

Wie groß ist n für die (n,n)-Koeffizientenmatrix?: 3

Werte der Koeffizientenmatrix eingeben:
Wert von a[0][0] eingeben: 1
Wert von a[0][1] eingeben: 1
Wert von a[0][2] eingeben: 1
Wert von a[1][0] eingeben: 2
Wert von a[1][1] eingeben: 2
Wert von a[1][2] eingeben: 2
Wert von a[2][0] eingeben: 1
Wert von a[2][1] eingeben: 1
Wert von a[2][2] eingeben: 1

Werte des Vektors b eingeben:
Wert von b[0] eingeben: 1
Wert von b[1] eingeben: 1
Wert von b[2] eingeben: 1

=========================================
|    x[0] |    x[1] |    x[2] |     b   |
=========================================
|       1 |       1 |       1 |       1 |
|       2 |       2 |       2 |       1 |
|       1 |       1 |       1 |       1 |
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|       1 |       1 |       1 |       1 |
|       0 |       0 |       0 |      -1 |
|       0 |       0 |       0 |       0 |
-----------------------------------------
=> 2. Zeile 0!=-1: also keine Lösung!

##########################################################################################
--------------------------------------
| Gauss-Algorithmus ohne Pivot-Suche |
--------------------------------------

Wie groß ist n für die (n,n)-Koeffizientenmatrix?: 3

Werte der Koeffizientenmatrix eingeben:
Wert von a[0][0] eingeben: 1
Wert von a[0][1] eingeben: 1
Wert von a[0][2] eingeben: 1
Wert von a[1][0] eingeben: 1
Wert von a[1][1] eingeben: 1
Wert von a[1][2] eingeben: 1
Wert von a[2][0] eingeben: 1
Wert von a[2][1] eingeben: 1
Wert von a[2][2] eingeben: 1

Werte des Vektors b eingeben:
Wert von b[0] eingeben: 1
Wert von b[1] eingeben: 1
Wert von b[2] eingeben: 1

=========================================
|    x[0] |    x[1] |    x[2] |     b   |
=========================================
|       1 |       1 |       1 |       1 |
|       1 |       1 |       1 |       1 |
|       1 |       1 |       1 |       1 |
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|       1 |       1 |       1 |       1 |
|       0 |       0 |       0 |       0 |
|       0 |       0 |       0 |       0 |
-----------------------------------------
=> 3. Zeile 0=0: also unendlich viele Lösungen!

##########################################################################################
--------------------------------------
| Gauss-Algorithmus ohne Pivot-Suche |
--------------------------------------

Wie groß ist n für die (n,n)-Koeffizientenmatrix?: 2

Werte der Koeffizientenmatrix eingeben:
Wert von a[0][0] eingeben: 5
Wert von a[0][1] eingeben: 3
Wert von a[1][0] eingeben: 4
Wert von a[1][1] eingeben: 6

Werte des Vektors b eingeben:
Wert von b[0] eingeben: 2
Wert von b[1] eingeben: 6

===============================
|    x[0] |    x[1] |     b   |
===============================
|       5 |       3 |       2 |
|       4 |       6 |       6 |
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|       5 |       3 |       2 |
|       0 |     3.6 |     4.4 |
===============================

Lösungen nach dem Gaussschen Algorithmus:
=> x[1] = 1.22222
=> x[0] = -0.333333									*/