/*　ガウス・ザイデル法による連立一次方程式の解  */

/*  ------ 各段階の近似値も表示する場合 ------  */

#include<stdio.h>

#include<math.h>

#define  N   4   /* number of eqations */

#define  EPSILON  1.0E-4

double  a[N][N]={{ 5.00, 2.00, 1.00, 2.55},

                { 1.00, 4.00, 2.00,-3.10},

                { 2.00, 1.00, 4.00, 0.97},

                {-0.20, 0.80,-0.25,11.01}};

double  b[N]={ 17.00, 16.00, 11.05, -8.25};

int    k;

double  x[N];

FILE *optfile;

void main()

{

    int i,j;

    void gaus_sei();

    optfile=fopen("08ガウス・ザイデル法.dat","w");

        fprintf(optfile,"====== Coefficient matrix and constant vector ====== \n\n");

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

    {

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

            fprintf(optfile,"%8.4f  ",a[i][j]);

        fprintf(optfile,"   %8.4f \n",b[i]);

    }

    fprintf(optfile,"\n==== Iteration of root (limit of error = %5.4f)====",EPSILON);

    fprintf(optfile,"\n\n");

    gaus_sei();

        fclose(optfile);

}

void gaus_sei()

{

    int i,j;

    double nk,sk,xx;   /* nk = error of n-th iteration */

    k=0;

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

        x[i]=0;       /* initial value (ｼｮｷﾁ) of root */

    fprintf(optfile," step");

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

        fprintf(optfile,"    x[%d]    ",j);

    fprintf(optfile,"\n");

    do

    {

        k=k+1;

        nk = 0.0;

        fprintf(optfile,"  %2d",k);

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

        {

           sk = 0.0;

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

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

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

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

           sk = b[i] - sk;

           xx = x[i];

           x[i] = sk / a[i][i];

           fprintf(optfile," %10.5f ",x[i]);

           nk = nk + fabs(xx - x[i]);

        }

        fprintf(optfile,"\n");

    }

    while(nk>EPSILON);

}