/*
 *  GauSed.c: GaussSeidelMethod ガウス・ザイデル法
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#ifndef max
#define max(a,b) ((a)>(b)?(a):(b))
#endif
#define CLINE printf(__FILE__":%d\n",__LINE__)

#define ND 10
#define IMAX 1000

int GauSed(double a[ND][ND], double b[ND], double x[ND], int n)
/*
  Gauss-Seidel method of solution of linear simultaneous equations
                   A * x = b
  Usage: GauSed(a, b, x, n);

    Arguments a ... [n*n] matrix
              b ... constant vector
              x ... solution vector
              n ... number of unknowns
*/
{
    double epsa=1.0e-70,epsr=1.0e-10;
    int i,j,itr,itrmx;
    double anorm,amax,xnew,xold,err,dx;

//  Check if this problem will converge. (収束予想)
//  itrmx is maximum iteration counter.
    anorm = 0.0;
    for (i=1;i<=n;i++)
    {
        amax = 0.0;
        for (j=1;j<=n;j++)
            if (i != j) amax = amax+fabs(a[i][j]);
        anorm = max(anorm, amax/fabs(a[i][i]));
    }
    if (anorm > 1.0)
    {
        putchar('\n');
        printf(" *** Gauss-Seidel method will not converge.\n");
        printf(" Execution stop on subroutine \"GauSed\".\n");
        printf(" Guessed matrix norm = %e\n",anorm);
        return 1;
    }
    if (anorm < 1) itrmx = (int)(log(epsr)/log(anorm));
    else           itrmx = IMAX;
    putchar('\n');
    printf(" Guessed matrix norm = %e\n",anorm);
    printf(" Maximum iteration counter = %d\n",itrmx);
    putchar('\n');

//  iteration (反復)
    itr = 0;
    do {
        itr++;
        err = 0.0;
        for (i=1;i<=n;i++)
        {
            xnew = 0.0;
            for (j=1;j<=i-1;j++) xnew += a[i][j]*x[j];
            for (j=i+1;j<=n;j++) xnew += a[i][j]*x[j];
            xnew = (b[i]-xnew)/a[i][i];
            xold = x[i];
            x[i] = xnew;
            dx = fabs(xnew-xold);
            err = max(err,dx/(epsa+epsr*(fabs(xnew)+fabs(xold))));
        }
        putchar('\n');
        printf("iteration =%d\n",itr);for (i=1;i<=n;i++) printf(" %e\n",x[i]);
    } while ((err >= 1.0) && (itr < itrmx));
    if (itr > itrmx)
    {
        printf(" *** Error *** Iteration does not converge.",
               " Execution terminated.\n");
        printf(" err = %e\n",err);
    }
    putchar('\n');printf(" number of iteration = %d\n",itr);
    return 0;
}

/*    ------------------    Main    ------------------    */
int main()
{
    int n,i,j;
    double a[ND][ND],b[ND],x[ND];
    char s[128],*s1;

    printf(" n = "); gets(s); sscanf(s,"%d",&n);
    if (n >= ND)
    {
      printf(" *** Input error *** n must be less than %2d.\n",ND);
      return 1;
    }
    for (i=1;i<=n;i++)
    {
        printf(" A[%2d][j] (j=1,%2d), B[%2d] = ",i,n,i);
        gets(s);
        s1=strtok(s," ");
        for (j=1;j<=n;j++) {
            sscanf(s1,"%lf",&a[i][j]);
            s1=strtok(NULL," ");
        }
        sscanf(s1,"%lf",&b[i]);
    }
    putchar('\n');

    printf("n = %d\n",n);

    printf("A =\n");
    for (i=1;i<=n;i++)
    {
        for (j=1;j<=n;j++) printf(" %e",a[i][j]);
        putchar('\n');
    }
    printf("b =\n");
    for (i=1;i<=n;i++) printf(" %e\n",b[i]);
    for (i=1;i<=n;i++) x[i] = 0;
// ----------------------
    GauSed(a, b, x, n);
// ----------------------
    putchar('\n');
    printf("Solution =\n");
    for (i=1;i<=n;i++) printf(" %e\n",x[i]);
    return 0;
}
/*   ---------   実行例: 教科書 p.97 例1 (4.60)式   ---------   */
/*
Z:\src\C>GauSed
 n = 2
 A[ 1,j] (j=1, 2), B[ 1] = 2 1 4
 A[ 2,j] (j=1, 2), B[ 2] = 1 2 5

n = 2
A =
 2.000000e+00 1.000000e+00
 1.000000e+00 2.000000e+00
b =
 4.000000e+00
 5.000000e+00

 Guessed matrix norm =  5.000000e-01
 Maximum iteration counter = 33

 number of iteration = 18

Solution =
 1.000000e+00
 2.000000e+00
*/