poisson_o.c

/*
 *  poisson_o:
 *
 *   solves the 2d poisson equation with simple Jacobi iterations 
 *   uses overlapping communication and computation
 *   source: simple quadrupole
 *
 *   usage:
 *              poisson N M
 *   
 *     N, M:  dimensions of global array A(N, M)
 *
 */

#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 and perform one Jacobi iteration */
    exchange_and_sweep1d(a, f, b, oned, grid1d);
    
    /* repeat with roles of a und b changed */
    error_local = exchange_and_sweep1d(b, f, a, oned, grid1d);

    /* 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);
}