/* * poisson.c * * solves the 2d poisson equation with simple Jacobi iterations * source: simple quadrupole * * usage: * poisson N M [IT] * * N, M: dimensions of global array A(N, M) * IT: max. number of iterations * */ #include "array.h" #include "poisson.h" void main(int argc, char *argv[]) { Decomposition oned; Grid grid1d; dmatrix a, b, f; int no_it; /* number of iterations */ int it; double error_local; double error_global; /* start tasks and broadcast infos */ no_it = startup(argc, argv, &oned, &grid1d); if (no_it < 1) { no_it = IT_MAX; } /* initialize the right hand side f and the initial guess a */ oned_init(&a, &b, &f, oned); /* do the iterations */ for (it=0; it<no_it; it++) { /* get the ghost points */ exchange(a, oned, grid1d); /* perform one Jacobi iteration */ sweep1d(a, f, b, oned); /* repeat with roles of a und b changed */ exchange(b, oned, grid1d); error_local = sweep1d(b, f, a, oned); /* check for convergence */ MPI_Allreduce(&error_local, &error_global, 1, MPI_DOUBLE, MPI_SUM, grid1d.comm); error_global /= oned.gx * oned.gy; if ((grid1d.me == 0) && (it % 100 == 0)) { printf("Iteration %3d: Difference is %10.6lf\n", it, error_global); } if (error_global < EPS) { break; } } if (grid1d.me == 0) { printf("\nGlobal error after %3d iterations: %10.6lf\n\n", it, error_global); } dump(a, oned, grid1d); finalize(a, b, f, &grid1d); }