exch-swp.c

/*
 * exchange-sweep.c
 *
 *   routine for one exchange-sweep-loop with overlapping of
 *      communication and computation
 *   replaces one consecutive call of exchange and sweep 
 */

#include "array.h"
#include "poisson.h"

double exchange_and_sweep1d(dmatrix a, dmatrix f, dmatrix b,
			    Decomposition oned, Grid grid) {
  int             lx, ly;
  MPI_Status      recv_status, send_status[2];
  MPI_Request     recv_req[2], send_req[2];
  int             i, j;
  int             index, it;
  double          error = 0.0;
  double          resid = 0.0;

  /*
   *   start all receives and sends  
   */
  lx  =  oned.lx;
  ly  =  oned.ly;
  
  /* extra boundary points in y are not exchanged, they are constant anyway */
  MPI_Irecv(&INDEX(a, 0, 1), ly, MPI_DOUBLE, grid.down,  1,
            grid.comm, &recv_req[0]);
  MPI_Irecv(&INDEX(a, lx+1, 1), ly, MPI_DOUBLE, grid.up, 0, 
            grid.comm, &recv_req[1]);
  MPI_Isend(&INDEX(a, 1, 1), ly, MPI_DOUBLE, grid.down, 0, 
	    grid.comm, &send_req[0]);
  MPI_Isend(&INDEX(a, lx, 1), ly, MPI_DOUBLE, grid.up,  1,
	    grid.comm, &send_req[1]);
  /*
   *   compute inner points
   */
  for (i = 2; i <= lx-1; i++) {
    for (j = 1; j <= ly; j++) {
      resid = INDEX(a, i-1, j) + INDEX(a, i, j+1) 
   	      + INDEX(a, i+1, j) + INDEX(a, i, j-1) - 4*INDEX(a, i, j)
	      - INDEX(f, i, j);
      error += fabs(resid);
      INDEX(b, i, j) = INDEX(a, i, j) + 0.25*resid;
    }
  }
  /*
   *  wait for the receives and compute corresponding edge points
   */ 
  for (it=0; it<=1; it++) {
    MPI_Waitany(2, recv_req, &index, &recv_status);
    
    /* which receive returned ? */
    if (index == 0) {
      i = 1;
    } else if (index == 1) {
      i = lx;
    } else {     /* index = MPI_UNDEFINED */
      printf("MPI_Waitany returned MPI_UNDEFINED!\n");
      MPI_Abort(MPI_COMM_WORLD, -1);
    }

    /* compute corresponding row */
    for (j = 1; j <= ly; j++) {
      resid = INDEX(a, i-1, j) + INDEX(a, i, j+1) 
  	      + INDEX(a, i+1, j) + INDEX(a, i, j-1) - 4*INDEX(a, i, j)
              - INDEX(f, i, j);
      error += fabs(resid);
      INDEX(b, i, j) = INDEX(a, i, j) + 0.25*resid;
    }
  }

  /*
   *   finally finish the sends
   */
  MPI_Waitall(2, send_req, send_status);
  
  return(error);
}