/*
 *  Jacobi.c: Jacobian Method ƒ„ƒRƒr–@
 */
#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 Jacobi(double a[][ND], double b[], double x[], double x0[], int n)
/*
  Jacobian method of solution of linear simultaneous equations
                   A * x = b
  Usage: Jacobi(a, b, x, x0, n);

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

//  Check if this problem will converge. (Žû‘©—\‘z)
//  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(" *** Jacobian method will not converge.\n");
        printf(" Execution stop on subroutine \"Jacobi\".\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++) x0[i] = x[i];
        for (i=1;i<=n;i++)
        {
            xnew = 0.0;
            for (j=1;j<=i-1;j++) xnew += a[i][j]*x0[j];
            for (j=i+1;j<=n;j++) xnew += a[i][j]*x0[j];
            xnew = (b[i]-xnew)/a[i][i];
            x[i] = xnew;
            dx = fabs(xnew-x0[i]);
            err = max(err,dx/(epsa+epsr*(fabs(xnew)+fabs(x0[i]))));
        }
        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);
    }
    printf(" number of iteration = %d\n",itr);
    return 0;
}

/*    ------------------    Main    ------------------    */
int main()
{
    int n,i,j;
    double a[ND][ND],b[ND],x[ND],x0[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;
// ------------------------
    Jacobi(a, b, x, x0, n);
// ------------------------
    putchar('\n');
    printf("Solution =\n");
    for (i=1;i<=n;i++) printf(" %e\n",x[i]);
    return 0;
}
/*   ---------   ŽÀ�s—á: ‹³‰È�‘ p.91 (4.45)Ž®   ---------   */
/*
Z:\src\C>Jacobi
 n = 2
 A[ 1][j] (j=1, 2), B[ 1] = 4 1 8
 A[ 2][j] (j=1, 2), B[ 2] = 1 2 5.5

n = 2
A =
 4.000000e+00 1.000000e+00
 1.000000e+00 2.000000e+00
b =
 8.000000e+00
 5.500000e+00

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

 number of iteration = 23

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