Jacobi C Code

#include 
int c100 = 100;
double cb3 = 1.5;
double cb7 = .8;
double cb8 = 1e-5;
int c20 = 20;

int  main()
{
    double uold[10000], f[10000], u[10000];
    double dx, dy;

/* *********************************************************** */
/* program to solve a finite difference */
/* discretization of Helmholtz equation : */
/* (d2/dx2)u + (d2/dy2)u - alpha u = f */
/* using Jacobi iterative method. */

/* Directives are used in this code to achieve paralleism. */
/* All do loops are parallized with 'even' scheduling to */
/* maximize performance. */

/* Input :  n - grid dimension in x direction */
/*          m - grid dimension in y direction */
/*          alpha - Helmholtz constant (always greater than 0.0) */
/*          tol   - error tolerance for iterative solver */
/*          relax - Successice over relaxation parameter */
/*          mits  - Maximum iterations for iterative solver */

/* On output */
/*       : u(n,m) - Dependent variable (solutions) */
/*       : f(n,m) - Right hand side function */
/* ************************************************************ */
/* ************************************************************* */
/* Working varaibles/arrays */
/*     dx  - grid spacing in x direction */
/*     dy  - grid spacing in y direction */
/* ************************************************************ */

/* Initialize data */

initialize_(&c100, &c100, &cb3, &dx, &dy, u, f);

/* Solve Helmholtz equation */

jacobi_(&c100, &c100, &dx, &dy, &cb3, &cb7, u, f, &cb8, &c20, uold);

/* Check error between exact solution */

error_check_(&c100, &c100, &cb3, &dx, &dy, u, f);
return (0);
} 

int initialize_(n, m, alpha, dx, dy, u, f)
int *n, *m;
double *alpha, *dx, *dy, *u, *f;
{
    int u_dim1, u_offset, f_dim1, f_offset, i1, i2;

    int ii, j, xx, yy;

/* ***************************************************** */
/* Initializes data */
/* Assumes exact solution is u(x,y) = (1-x^2)*(1-y^2) */

/* ***************************************************** */
    /* Parameter adjustments */
    f_dim1 = *n;
    f_offset = f_dim1 + 1;
    f -= f_offset;
    u_dim1 = *n;
    u_offset = u_dim1 + 1;
    u -= u_offset;

    /* Function Body */
    *dx = (float)2. / (*n - 1);
    *dy = (float)2. / (*m - 1);
/* Initilize initial condition and RHS */
    i1 = *m;
    for (j = 1; j <= i1; ++j) {
	i2 = *n;
	for (ii = 1; ii <= i2; ++ii) {
	    xx = (int) (*dx * (double) (ii - 1) - (float)1.);

/* -1 < x < 1 */

	    yy = (int) (*dy * (double) (j - 1) - (float)1.);

/* -1 < y < 1 */

	    u[ii + j * u_dim1] = (float)0.;
	    f[ii + j * f_dim1] = -(*alpha) * ((float)1. - xx * xx) * ((float)
		    1. - yy * yy) - ((float)1. - xx * xx) * (float)2. - ((
		    float)1. - yy * yy) * (float)2.;
	}
    }
    return 0;
} /* initialize_ */

int jacobi_(n, m, dx, dy, alpha, omega, u, f, tol, maxit, uold)
int *n, *m;
double *dx, *dy, *alpha, *omega, *u, *f, *tol;
int *maxit;
double *uold;
{
    /* System generated locals */
int f_dim1, f_offset, u_dim1, u_offset, uold_dim1, uold_offset, i1, i2;

    /* Builtin functions */
    double sqrt();
    int s_wsle(), do_lio(), e_wsle();

    /* Local variables */
    double b;
    int ii, j, k;
    double resid, error, ax, ay;

/* ***************************************************************** */
/* Subroutine HelmholtzJ */
/* Solves poisson equation on rectangular grid assuming : */
/* (1) Uniform discretization in each direction, and */
/* (2) Dirichlect boundary conditions */

/* Jacobi method is used in this routine */

/* Input : n,m   Number of grid points in the X/Y directions */
/*         dx,dy Grid spacing in the X/Y directions */
/*         alpha Helmholtz eqn. coefficient */
/*         omega Relaxation factor */
/*         f(n,m) Right hand side function */
/*         u(n,m) Dependent variable/Solution */
/*         tol    Tolerance for iterative solver */
/*         maxit  Maximum number of iterations */

/* Output : u(n,m) - Solution */
/* **************************************************************** */

/* Local variables */


/* Initialize coefficients */
    /* Parameter adjustments */
    uold_dim1 = *n;
    uold_offset = uold_dim1 + 1;
    uold -= uold_offset;
    f_dim1 = *n;
    f_offset = f_dim1 + 1;
    f -= f_offset;
    u_dim1 = *n;
    u_offset = u_dim1 + 1;
    u -= u_offset;

    /* Function Body */
    ax = (float)1. / (*dx * *dx);

/* X-direction coef */

    ay = (float)1. / (*dy * *dy);

/* Y-direction coef */

    b = (float)-2. / (*dx * *dx) - (float)2. / (*dy * *dy) - *alpha;

/* Central coeff */

    error = *tol * (float)10.;
    k = 1;
    while(k <= *maxit && error > *tol) {
	error = (float)0.;
/* Copy new solution into old */
	i1 = *m;
	for (j = 1; j <= i1; ++j) {
	    i2 = *n;
	    for (ii = 1; ii <= i2; ++ii) {
		uold[ii + j * uold_dim1] = u[ii + j * u_dim1];
	    }
	}
/* Compute stencil, residual, & update */
	i1 = *m - 1;
	for (j = 2; j <= i1; ++j) {
	    i2 = *n - 1;
	    for (ii = 2; ii <= i2; ++ii) {
/*     Evaluate residual */
		resid = (ax * (uold[ii - 1 + j * uold_dim1] + uold[ii + 1 + 
			j * uold_dim1]) + ay * (uold[ii + (j - 1) * 
			uold_dim1] + uold[ii + (j + 1) * uold_dim1]) + b * 
			uold[ii + j * uold_dim1] - f[ii + j * f_dim1]) / b;
/* Update solution */
		u[ii + j * u_dim1] = uold[ii + j * uold_dim1] - *omega * 
			resid;
/* Accumulate residual error */
		error += resid * resid;
	    }
	}
/* Error check */
	++k;
	error = sqrt(error) / (double) (*n * *m);

    }

printf("Total Number of Iterations %d\n" , k);
printf("Residual                   %g\n" ,error);
    return 0;
} /* jacobi_ */

int error_check_(n, m, alpha, dx, dy, u, f)
int *n, *m;
double *alpha, *dx, *dy, *u, *f;
{
    /* System generated locals */
    int u_dim1, u_offset, f_dim1, f_offset, i1, i2;

    /* Builtin functions */
    double sqrt();

    /* Local variables */
    double temp;
    int ii, j;
    double error, xx, yy;


/* *********************************************************** */
/* Checks error between numerical and exact solution */

/* *********************************************************** */
    /* Parameter adjustments */
    f_dim1 = *n;
    f_offset = f_dim1 + 1;
    f -= f_offset;
    u_dim1 = *n;
    u_offset = u_dim1 + 1;
    u -= u_offset;

    /* Function Body */
    *dx = (float)2. / (*n - 1);
    *dy = (float)2. / (*m - 1);
    error = (float)0.;
    i1 = *m;
    for (j = 1; j <= i1; ++j) {
	i2 = *n;
	for (ii = 1; ii <= i2; ++ii) {
	    xx = *dx * (double) (ii - 1) - 1.;
	    yy = *dy * (double) (j - 1) - 1.;
	    temp = u[ii + j * u_dim1] - ((float)1. - xx * xx) * ((float)1. - 
		    yy * yy);
	    error += temp * temp;
	}
    }
    error = sqrt(error) / (double) (*n * *m);
    printf("Solution Error : %g\n" , error);
    printf("i =    %d\n", ii);
    printf("j =    %d\n", j);
    return 0;
} /* error_check__ */