# 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) */

/* ***************************************************** */
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 */
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 */

/* *********************************************************** */
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__ */