poisson.c

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