/*  ガウスの消去法による連立方程式の解　*/

#include<stdio.h>

#include<math.h>

#define  N      3  /* number of equations */

#define  EPSILON  1.0E-5

#define  TRUE   1

#define  FALSE  0

float  a[N][N]={{ 2.0, 5.0,-3.0},

                { 4.0, 5.0,-3.0},

                {-1.0,-5.0, 4.0}};

float  x[N],b[N]={ 13.0, 19.0,-9.0};

int    sw;

void main()

{

    void gauss(int *sw);

    int i,j;

    printf("\n\n ===== Coefficient matrix and constant vector ===== \n\n");

    for(i=0;i<N;i++)

    {

        for(j=0;j<N;j++)

            printf("  %10.6f",a[i][j]);

        printf("     %10.6f\n",b[i]);

    }

    sw=TRUE;

    gauss(&sw);

    printf("\n\n ==== Coefficient matrix and constant vector \n");

    printf("                   by Gauss elimination method ====\n\n");

    for(i=0;i<N;i++)

    {

        for(j=0;j<N;j++)

            printf("  %10.6f",a[i][j]);

        printf("     %10.6f\n",b[i]);

    }

    printf("\n\n ========== Root of simultaneous eqations ========= \n\n");

    if(sw=TRUE)

        for(i=0;i<N;i++)

            printf("  x[%d] = %12.6f \n",i+1,x[i]);

    else

        printf("        can't be find.\n");

    printf("\n");

}

void swap(float *wk1,float *wk2)

{

    float w;

    w = *wk1;

    *wk1 = *wk2;

    *wk2 = w;

}

void gauss(int *sw)

{

    int i,j,k,m,ik;

    float aik,akk,sa;

    for(k=0;k<N;k++)

    {

        akk = a[k][k];

        if(fabs(akk)<=EPSILON)

        {

            ik = k+1;

            while((ik<N)&&(fabs(a[ik][k])<EPSILON))

                ik++;

            if(ik<N)

            {

                for(m=k;m<N;m++)

                    swap(&a[k][m],&a[ik][m]);

                swap(&b[k],&b[ik]);

                akk = a[k][k];

            }

            else

            {

                printf("    pivot = 0  ! \n");

                *sw=FALSE;

                goto end;

            }

        }

        for(j=k;j<N;j++)

            a[k][j]=a[k][j]/akk;

        b[k]=b[k]/akk;

        if(k<N-1)

        {

            for(i=k+1;i<N;i++)

            {

                aik=a[i][k];

                for(j=k;j<N;j++)

                    a[i][j]=a[i][j]-a[k][j]*aik;

                b[i]=b[i]-b[k]*aik;

            }

        }

    }

    x[N-1]=b[N-1];

    for(i=N-2;i>=0;i--)

    {

        sa=0;

        for(j=i+1;j<N;j++)

            sa=sa+a[i][j]*x[j];

        x[i]=b[i]-sa;

    }

    end:;

}