/*
 *  Gauss.c GaussElimMethod ガウスの消去法
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define ND 10
#define MD 20
#define max(a,b) ((a)>(b)?(a):(b))

int Gauss(double a[ND][MD], double *det, int n, int m)
/*
  Gaussian eliminate method of solution of linear
  simultaneous equations with partial pivotting for
                   A * x = b
  Usage: Gauss(a, &det, n, m);

    Arguments  a ..... A[n][n+m] matrix
               det ... determinant of A[n][m]
               n ..... number of unknowns
               m ..... number of problems
    Note: b[i][j] is stored in a[i][n+j]  (1<=j<=m)
          solution x[i][j] is stored in a[i][n+j]  (1<=j<=m)
*/
{
    int i, j, k, L;
    double eps, w, p, aik;

  /* initialization (初期化) */
    *det = 1.0; eps = 0.0;
    for(k=1;k<=n;k++)
        for(j=1;j<=n;j++) eps = max(eps,fabs(a[k][j]));
    eps = eps*1.0e-10;

/* forward elimination (前進消去) */
    for(k=1;k<=n;k++)
    {
    /* find pivot (ピボット選択) */
    /* - partial pivotting (部分ピボット選択法) */
        L = k;
        for(j=k+1;j<=n;j++)
            if(fabs(a[k][k]) < fabs(a[j][k])) L = j;
    /* exchange rows if necessary (行の交換) */
        if(L != k)
        {
            *det = -*det;   /* 行列式の符号反転 */
            for(j=k;j<=n+m;j++)
            {
                w = a[k][j]; a[k][j] = a[L][j]; a[L][j] = w;
            }
        }
    /* end of pivotting (ピボット選択終了) */

        p = a[k][k];

    /* determinant (行列式の計算) */
        *det = *det*p;

    /* return if matrix is singular (特異行列の処理) */
        if(fabs(*det) < eps)
        {
            perror(" *** error *** matrix is singular.");
            exit(1);
        }

        for(j=k+1;j<=n+m;j++) a[k][j] = a[k][j]/p;
        for(i=k+1;i<=n;i++)
        {
            aik = a[i][k];
            for (j=k+1;j<=n+m;j++) a[i][j] = a[i][j]-aik*a[k][j];
        }
    }

/* backword substitution (後退代入) */
    for(k=n-1;k>=1;k--)
        for(i=k+1;i<=n;i++)
            for(j=n+1;j<=n+m;j++)
                a[k][j] = a[k][j]-a[k][i]*a[i][j];
    return EXIT_SUCCESS;
}

/*    ------------    Main    ------------    */
int main()
{
    int m,n,i,j;
    double det;
    double a[ND][MD];
    char s[128],*s1;

    printf(" n = "); gets(s); sscanf(s,"%d",&n);
    printf(" m = "); gets(s); sscanf(s,"%d",&m);

    if((n > ND) | (n+m > MD))
    {
        printf(" *** input error *** n and n+m must be less",
             " than %2d and %2d, respectively.",ND,MD);
        exit(EXIT_FAILURE);
    }

    for(i=1;i<=n;i++)
    {
        printf(" A[%2d][j] (j=1,%2d) = ",i,n+m);
        gets(s);
        s1=strtok(s," ");
        for(j=1;j<=n+m;j++)
        {
            sscanf(s1,"%lf",&a[i][j]);
            s1=strtok(NULL," ");
        }
    }
    putchar('\n');

    printf("A =\n");
    for(i=1;i<=n;i++)
    {
        for(j=1;j<=n;j++) printf(" %e",a[i][j]);
        putchar('\n');
    }
    putchar('\n');

    printf("b =\n");
    for(i=1;i<=n;i++)
    {
        for(j=n+1;j<=n+m;j++) printf(" %e",a[i][j]);
        putchar('\n');
    }
    putchar('\n');
/* --------------------- */
    Gauss(a, &det, n, m);
/* --------------------- */
    printf("x =\n");
    for(i=1;i<=n;i++)
    {
        for(j=n+1;j<=n+m;j++) printf(" %e",a[i][j]);
        putchar('\n');
    }
    putchar('\n');
    printf("determinant = %e\n", det);
    return EXIT_SUCCESS;
}

/*   ---------   実行例: 教科書 p.54 (4.1)式   ---------   */
/*
Z:\src\C>Gauss
 n = 3
 m = 1
 A[ 1][j] (j=1, 4) = 3 2 1 10
 A[ 2][j] (j=1, 4) = 2 4 -1 7
 A[ 3][j] (j=1, 4) = 1 1 5 18

A =
 3.000000e+00 2.000000e+00 1.000000e+00
 2.000000e+00 4.000000e+00 -1.000000e+00
 1.000000e+00 1.000000e+00 5.000000e+00

b =
 1.000000e+01
 7.000000e+00
 1.800000e+01

x =
 1.000000e+00
 2.000000e+00
 3.000000e+00

determinant = 3.900000e+01
*/