/*  Jacobi method (roots of simultaneous equations) */

/*  Gauss Seidel method ÄÉ Ë¶¸  */

#include<stdio.h>

#include<math.h>

#define  N  3

#define  EPSILON 1.0E-4

double a[N][N]={{ 5.0, 2.0, 1.0 },

               { 1.0, 4.0, 2.0 },

               { 2.0, 1.0, 4.0 }};

double b[N]={ 17.0, 16.0, 11.0 };

double xJ[N],xG[N];

int k=0;  /* k-th stage of iteration */

double xJ[N],xG[N];  /* roots(Jacobi,Gauss Seidel) */

int swJ=0,swG=0,swpJ=0,swpG=0;

FILE *optfile;

void main()

{

    int i;   

        void Jacobi(void);

    void GaussSeidel(void);

        void opt(void);

    void opt_roots(void);

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

        xJ[i]=xG[i]=0.0;

        optfile=fopen("07ƒ„ƒRƒr–@.dat","w");

    opt();

    opt_roots();

    while(swJ==0 || swG==0)

    {

        k=k+1;

        if(swJ==0)

            Jacobi();

        if(swG==0)

            GaussSeidel();

        opt_roots();

    }

        fclose(optfile);

}

void Jacobi()

{

    int i,j;

    double e,s,x[N];

    e=0.0;

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

    {

        s=b[i];

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

        {

            if(i!=j)

                s=s-a[i][j]*xJ[j];

        }

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

        e=e+fabs(x[i]-xJ[i]);

    }

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

        xJ[i]=x[i];

    if(e<=EPSILON)

        swJ=1;

}

void GaussSeidel()

{

    int i,j;

    double e,s,x;

    e=0.0;

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

    {

        x=xG[i];

        s=b[i];

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

            if(i!=j)

                s=s-a[i][j]*xG[j];

        xG[i]=s/a[i][i];

        e=e+fabs(x-xG[i]);

    }

    if(e<=EPSILON)

        swG=1;

}

void opt()

{

    int i,j;

    fprintf(optfile,"------ Coeficient matrix & constant vector ------\n\n");

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

    {

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

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

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

    }

    fprintf(optfile,"\n----------------------- Roots of k-th stage ");

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

    fprintf(optfile,"stage   ****** Jacobi method ******");

    fprintf(optfile,"     *** Gauss Seidel method ***\n\n");

    fprintf(optfile,"  #        x1        x2        x3");

    fprintf(optfile,"          x1        x2        x3 \n");

}

void opt_roots()

{

    int i;

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

    swpJ=swpJ+swJ;swpG=swpG+swG;

    if(swpJ<=1 && swpG<=1)

    {

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

            fprintf(optfile,"  %8.5f",xJ[i]);

        fprintf(optfile,"  ");

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

            fprintf(optfile,"  %8.5f",xG[i]);

    }

    if(swpJ<=1 && swpG>1)

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

            fprintf(optfile,"  %8.5f",xJ[i]);

}