/*
 *  lineq.c : Lösung eines n x n linearen Gleichungssystems
 *            A*x=b mit dem Gaußschen Eliminationsverfahren
 *            ohne Pivotsuche.
 *
 *  a[1..n][1..n] ... n x n Koeffizientenmatrix
 *  b[1..n] ......... rechte Seite des Systems
 *  x[1..n] ......... Lösungsvektor des Systems
 */

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

#define  EPS  1.0e-8

/*
 *  Gaußsches Eliminationsverfahren ohne Pivotsuche
 */
int gauss( int n, double **a, double *b, double *x ) {

    double pivot, zmult, sum;
    int    i, j, k;

    /*  Elimination  */

    for (i = 1; i <= n - 1; i++) {      /* Pivotzeile                      */

        pivot = a[i][i];                /* Pivotelement                    */

        if (fabs(pivot) <= EPS) {
           return(1);
        }

        for (k = i + 1; k <= n; k++) {  /* Eliminationsschritt für Zeile k */

            zmult = a[k][i]/pivot;

            for (j = i + 1; j <= n; j++) {    /* neues Element in Spalte j */
                a[k][j] = a[k][j] - zmult*a[i][j];
            }

            b[k] = b[k] - zmult*b[i];
        }
    }

    /*  Rückwärtseinsetzen  */

    if (fabs(a[n][n]) <= EPS) {     /* alle anderen a[i][i] schon getestet */
       return(1);
    }

    x[n] = b[n]/a[n][n];

    for (i = n - 1; i >= 1; i--) {

        sum = b[i];

        for (j = i + 1; j <= n; j++) {
            sum = sum - a[i][j]*x[j];
        }

        x[i] = sum/a[i][i];
    }

    return(0);
}


int main( void ) {

    double **a;
    double *b, *x;
    int    i, j, n;

    printf("n: ");
    scanf("%d", &n);

    a = (double **) calloc( n+1, sizeof(double *) );
    for (i = 0; i <= n; i++) {
        a[i] = (double *) calloc( n+1, sizeof(double) );
    }

    b = (double *) calloc( n+1, sizeof(double) );
    x = (double *) calloc( n+1, sizeof(double) );

    /*  Koeffizientenmatrix a und rechte Seite b einlesen  */
    for (i = 1; i <= n; i++) {
        for (j = 1; j <= n; j++) {
            scanf("%lf", &a[i][j]);
        }
        scanf("%lf", &b[i]);
    }

    /*  Lösungsvektor ausgeben  */
    if (gauss(n, a, b, x) == 0) {
       printf("\nx = ( ");
       for (i = 1; i <= n - 1; i++) {
           printf("%f, ", x[i]);
       }
       printf("%f )\n", x[n]);
    } else {
       fprintf(stderr, "lineq: Gauss-Elimination nicht durchfuehrbar\n");
       exit(1);
    }

    return(0);
}