/*
 *  lineqp.c: Lösung eines n x n linearen Gleichungssystems
 *            A*x=b mit dem Gaußschen Eliminationsverfahren
 *            mit Spaltenpivotsuche.
 *
 *  a[1..n][1..n] .. n x n Koeffizientenmatrix
 *  b[1..n] ........ rechte Seite des Systems
 *  x[1..n] ........ Lösungsvektor des Systems
 *  p[1..n] ........ Permutationsvektor zur Buchführung der
 *                   Zeilenvertauschungen
 */


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

#define  EPS  1.0e-8


/*
 *  Gaußsches Eliminationsverfahren mit Spaltenpivotsuche
 */
int gaussp( int n, double **a, double *b, double *x, int *p ) {

    double pivot, pmax, ptemp, zmult, sum;
    int    i, j, k, m, pi, pk, itemp;


    /*  Permutationsvektor initialisieren  */
    for (i = 1; i <= n; i++)
        p[i] = i;


    /*  Elimination  */

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

        /*  Pivotsuche  */

        m = i;
        pmax = 0.0;
        for (j = i; j <= n; j++) {    /* suche max|a[j][i]| unterhalb     */
            ptemp = fabs(a[p[j]][i]); /* der Hauptdiagonale in Spalte i   */
            if (ptemp > pmax) {
               pmax = ptemp;
               m = j;                 /* Zeilenindex für neue Pivotzeile  */
            }
        }

        if ( pmax <= EPS ) {          /* Matrix nahezu singulär           */
           return(1);
        }

        if ( m != i ) {               /* Zeile i mit Zeile m vertauschen  */
           itemp = p[i];
           p[i] = p[m];
           p[m] = itemp;
        }

        /*  Eliminationsschritte für verbleibende Zeilen k = i+1,...,n    */

        pi = p[i];
        pivot = a[pi][i];

        for (k = i + 1; k <= n; k++) {

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

            /*  Neue Elemente in den Spalten j der Zeile k  */
            for (j = i + 1; j <= n; j++)
                a[pk][j] = a[pk][j] - zmult*a[pi][j];

            b[pk] = b[pk] - zmult*b[pi];
        }
    }

    /*  Rückwärtseinsetzen  */

    if ( fabs(a[p[n]][n]) <= EPS ) {
       return(1);
    }

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

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

        pi = p[i];
        sum = b[pi];

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

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

    return(0);
}



int main( void ) {

    double **a;
    double *b, *x;
    int    *p;
    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) );
    p = (int *) calloc( n+1, sizeof(int) );

    /*  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 (gaussp(n, a, b, x, p) == 0) {
       printf("\nx = ( ");
       for (i = 1; i <= n - 1; i++)
           printf("%f, ", x[i]);
       printf("%f )\n", x[n]);
    } else {
       fprintf(stderr, "lineqp: Koeffizientenmatrix singulaer\n");
       exit(1);
    }

    return(0);
}