/*
 *  Gaujor.c Gauss-Jordan Method
 *  ガウス・ジョルダンの消去法(掃出法)
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define ND 10
#define MD 20
#ifndef max
#define max(a,b) ((a)>(b)?(a):(b))
#endif

int GauJor(double a[ND][MD], double *det, int n, int m)
/*
  Gauss-Jordan eliminate (sweep out) method of solution of
  linear simultaneous equations with partial pivotting for
                   A * x = b
  Usage: GauJor(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(EXIT_FAILURE);
        }

        for (j=k+1;j<=n+m;j++) a[k][j] /= p;
        for (i=1;i<=n+1;i++)
            if (i != k)
            {
                aik = a[i][k];
                for (j=k+1;j<=n+m;j++) a[i][j] -= aik*a[k][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');
// -----------------------
    GauJor(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.62 (4.19)式   ---------   */
/*
Z:\src\C>GauJor
 n = 3
 m = 1
 A[ 1][j] (j=1, 4) = 4 1 1 9
 A[ 2][j] (j=1, 4) = 1 3 1 10
 A[ 3][j] (j=1, 4) = 2 1 5 19

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

b =
 9.000000e+00
 1.000000e+01
 1.900000e+01

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

determinant = 4.800000e+01
*/