/*
 *  lineql.c : Lösung eines n x n linearen Gleichungssystems
 *             A*x=b mit der LAPACK-Routine dgesv().
 *
 *  a[1..n][1..n] ... n x n Koeffizientenmatrix
 *  b[1..n] ......... rechte Seite des Systems
 *
 *  Linux:
 *  cc -o lineql lineql.c -L/usr/lib/lapack -llapack -lblas -lg2c -lm
 */

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

#define   EPS  1.0e-8

extern void dgesv_( int *n, int *nrhs, double *at, int *lda, int *p,
                    double *b, int *ldb, int *info );


int main( void ) {

    double **a;
    double *b;
    double *at;
    int    *p;
    int    i, j, n;
    int    nrhs, lda, ldb, info;

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

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

    b = (double *) calloc( n, sizeof(double) );
    p = (int *) calloc( n, sizeof(int) );

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

    /*  In Fortran wird eine Matrix  s p a l t e n w e i s e    */
    /*  in einem eindimensionalen Feld gespeichert.             */
     
    at = (double *) malloc( (size_t) ((n * n) * sizeof(double)) );
    for (j = 0; j <= n - 1; j++) {
        for (i = 0; i <= n - 1; i++) {
            at[j*n + i] = a[i][j];
        }
    }

    nrhs = 1;               /* number of right hand sides */
    lda  = n;               /* leading dimension of a     */
    ldb  = n;               /* leading dimension of b     */

    (void) dgesv_( &n, &nrhs, at, &lda, p, b, &ldb, &info );

    /*  in b wird der Lösungsvektor zurückgegeben, in  at */
    /*  die Faktoren L und U der LU-Zerlegung in  a=P*L*U */

    if (info == 0) {
       /*  Matrix  a  numerisch singulär?                             */
       /*  dgesv() testet nur auf   e x a k t e   Singularität !      */
       for (i = 0; i <= n - 1; i++) {
           if (fabs(at[i*n + i]) <= EPS) {       /*  U[i,i] <= EPS ?  */
              fprintf(stderr, "lineql: Matrix (numerisch) singulaer\n");
              exit(1);
           }
       }
       /* Lösung ausgeben, wenn Matrix a  regulär */
       printf("\nx = ( ");
       for (i = 0; i <= n - 2; i++)
           printf("%f, ", b[i]);
       printf("%f )\n", b[n-1]);
    } else {
       fprintf(stderr, "lineql: info = %d\n", info);
       exit(1);
    }

    return(0);
}