/*
 *  LUdec.c -- LU Factorization Method (LU分解法)
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define ND 10
#define MD 10
#ifndef max
#define max(a,b) ((a)>(b)?(a):(b))
#endif

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

    Arguments  a ..... A[n][n] matrix
               b ..... constant vector
               x ..... solution vector
               det ... determinant of A[n][m]
               n ..... number of unknowns
               m ..... number of problems
               mm .... order of pivotting (並べ換えリスト)
*/
{
    int i, j, k, L, iw;
    double w, p, aik;

//  LU factorization (LU分解)

    for (i=1;i<=n;i++) mm[i] = i;
    *det = 1.0;

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

        p = a[k][k];

        for (j=k+1;j<=n;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;j++) a[i][j] = a[i][j]-aik*a[k][j];
        }
    }

//  determinant (行列式の計算)
    for (k=1;k<=n;k++) *det *= a[k][k];
    return EXIT_SUCCESS;
}

int Solve(double a[][ND], double b[][MD], double x[][MD], int n, int m, int mm[])
{
    int i, j, k;

    for (j=1; j<=m;j++)
    {
    //  forward substisution (前進代入)
        for (k=1;k<=n;k++)
        {
            x[k][j] = b[mm[k]][j]; // Sort in pivotting order
            for (i=1;i<=k-1;i++) x[k][j] = x[k][j]-a[k][i]*x[i][j];
            x[k][j] = x[k][j]/a[k][k];
        }
    //  backward substisution (後退代入)
        for (k=n;k>=1;k--)
            for (i= k+1;i<=n;i++) x[k][j] = x[k][j]-a[k][i]*x[i][j];
    }
    return EXIT_SUCCESS;
}

/*    ------------    Main    ------------    */
int main()
{
    int m,n,i,j,mm[ND];
    double det;
    double a[ND][ND],b[ND][MD],x[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) | (m >= MD))
    {
        printf(" *** input error *** n and 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), b[%2d][j] (j=1,%2d) = ",i,n,i,m);
        gets(s);
        s1=strtok(s," ");
        for (j=1;j<=n;j++)
        {
            sscanf(s1,"%lf",&a[i][j]);
            s1=strtok(NULL," ");
        }
        for (j=1;j<=m;j++)
        {
            sscanf(s1,"%lf",&b[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=1;j<=m;j++) printf(" %e",b[i][j]);
        putchar('\n');
    }
    putchar('\n');
// ---------------------
    LU(a, &det, n, mm);
// ---------------------
    printf("Factorized A =\n");
    for (i=1;i<=n;i++)
    {
        for (j=1;j<=n;j++) printf(" %e",a[i][j]);
        putchar('\n');
    }
    putchar('\n');

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

/*   ---------   実行例: 教科書 p.88 (4.42)式   ---------   */
/*
Z:\src\C>LUdec
 n = 3
 m = 1
 A[ 1,j] (j=1, 3), b[ 1,j] (j=1, 1) = 3 2 1 6
 A[ 2,j] (j=1, 3), b[ 2,j] (j=1, 1) = 2 4 3 10
 A[ 3,j] (j=1, 3), b[ 3,j] (j=1, 1) = 1 3 5 12.5

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

b =
 6.000000e+00
 1.000000e+01
 1.250000e+01

Factorized A =
 3.000000e+00 6.666667e-01 3.333333e-01
 2.000000e+00 2.666667e+00 8.750000e-01
 1.000000e+00 2.333333e+00 2.625000e+00

determinant = 2.100000e+01

x =
 1.000000e+00
 5.000000e-01
 2.000000e+00
*/