/*    !!! Hausaufgabe: Ergaenzen Sie die Ausgabeprozedur 'write_x"   !!!       */

/* Gauss-Elimination mit Gesamtpivot-Suche. Die Koeffizienten muessen in einer */
/* Datei abgelegt sein i. d. Reihenfolge: 			               */
/* a11 a12 ... a1N b1 a21 a22 ... a2N b2 ...... aNN bN                         */
/* Zwischen die Werte sind Blanks zu setzen; Zeilenumbrueche sind  erlaubt.    */
/* Das Programm erfragt den Dateinamen.                                        */

#include <stdio.h>
#include <math.h>
#define N 3                /* Dimension der Koeffizienten-Matrix    */
double a[N+1][N+2];        /* Erweiterte Koeffizienten-Matrix       */
double x[N+1];             /* Loesungsvektor                        */
void read_a(void);         /* Einlesen d. erweiterten Koeff.-Matrix */  
void write_x(void);        /* Loesungsvektor ausgeben               */
void eliminate(void);
void substitute(void);
void gesamtpivot(int i);    
void write_a(void);        /* Erw. Koeffizienten-Matrix ausgeben    */

main()
{
	read_a();
	write_a();
	eliminate();
	substitute();
	write_x();
}

void read_a(void)
/* Einlesen der erweiterten Koeff.-Matrix aus einer Datei */
{
int z, s;     /* z = Zeilenindex, s = Spaltenindex */
/* Zeile 0 kontrolliert die Spaltenvertauschungen */
/* Spalte 0 wird nicht verwendet */
FILE *data_in;
char datname[15];
do {printf("Geben Sie den Namen der Datei mit den Matrix-Koeff. ein : ");
	scanf("%[^\n]", datname);
	data_in = fopen(datname, "r");
	if (data_in == NULL) printf("Fehler!\n");
   } while (data_in == NULL);
for (s=1; s<=N; s++) a[0][s] = s; /* Verwaltung der Indices bei Gesamtpivotsuche */
for (z = 1; z<= N; z++)
	for (s = 1; s <= N+1; s++)
	fscanf(data_in, "%lf", &a[z][s]);
fclose(data_in);
}

void eliminate(void)
/* Koeff.-Matrix in 3-Ecksgestalt bringen */
{
int i, j, k;
double piv, q;
for (i=1; i<N; i++)
	{
	gesamtpivot(i);
	piv = a[i][i];
	for (j=i+1; j<=N; j++)
		{
		q=a[j][i]/piv;
		for (k=i; k<=N+1; k++)
		    a[j][k] -= a[i][k]*q;
		}
	}
}

void substitute(void)
/* 3-ecksmatrix nach x auflösen; Indexe sind nicht korrekt;*/
/* Korrektur erfolgt anhand Zeile 0 in Funktion writex */
{
int j, k;
double t;
for (j = N; j>=1; --j)
	{ t=0.0;
	for (k=j+1; k<=N; ++k)
	    t += a[j][k]*x[k];
	x[j] = (a[j][N+1]-t)/a[j][j];
write_a();
	}

}





/**********************************************************************/
void write_x(void)
/* Ausgabe des korrekten Lösungsvektors anhand Zeile 0 */
{
  
  int i, j;
  // i = Index der auszugebenden Variablen
  for(i=1; i<=N; i++)
    for(j=1; j<=N; j++)
      //Ueberpruefen der Zeile 0 von a, da diese die Indizes enthaelt
      if(static_cast<int>(a[0][j])==i){
	//Ausgeben von xi und dem, zum Index in a passenden,
        //Wert aus x
        printf(" x%d=%15.6f", i, x[j]);
	printf("\n");
      }  	
        
    

/* Die Ausgabe muss in der folgenden Form erscheinen:		      */	
/*	x1 =     						      */
/*      x2 =    						      */	
/*         .							      */
/*         .							      */	
/*	xN =     						      */
/* Das muss auch fuer N > 3 funktionieren!			      */
	
}	
/**********************************************************************/

void gesamtpivot(int i)
/* Pivotsuche in der Untermatrix ab Diagonalelement a[i][i] */
{
double t;
int z, s, maxz, maxsp;
maxz = maxsp = i;
for (z=i; z<=N; ++z)
	for (s=i; s<=N; ++s)
		if (fabs(a[z][s])>fabs(a[maxz][maxsp]))
			{maxz = z;
			 maxsp = s;
			}
/* Vertausche Zeile maxz mit Zeile i */
for (s=i; s<=N+1; ++s)
	{t=a[i][s];
	   a[i][s]=a[maxz][s];
	           a[maxz][s]=t;
	}
/* Vertausche Spalte maxsp mit Spalte i */
for (z=0; z<=N; ++z)
	{t=a[z][i];
	   a[z][i]=a[z][maxsp];
	           a[z][maxsp]=t;
	}
}

void write_a(void)
/* Ausgabe der gesamten erweiterten Koeff.-Matrix */
/* Funktion an beliebigen Stellen d. Progr. einfügbar */
{
int z, s;
int c;
for (s=1; s<=N; ++s)
	printf(" %7.0f   ", a[0][s]);
printf("\n");
for (z=1; z<=N; ++z)
	{
	for (s=1; s<=N+1; ++s)
	    printf(" %10.6f", a[z][s]);
	printf("\n");
	}
}