redblack2d.c

LPARX is a C++ class library for for implementing portable scientific applications on distributed me

//
// Perform Red-Black iterations for two dimensional Laplace's equation
// This is intended to be an example of how to use LPARX.
//

#include "redblack2d.h"

/*
  Each LPARX grid has two attributes:

	1. the codomain, the base type. e.g. double
	2. the rank, the number of dimensions

  LPARX codomains can be of any type.
  The current implementation supports ranks 1 through 4.

  The name of a Grid type is formed as follows

	GridX(T)

  where X is the rank (an int in the range 1..4) and T is a type.

  Macros for types Grid1, Grid2, Grid3, and Grid4 are defined in
  the standard header file "lparx.h", which must be included
  in every LPARX module (this include is found in "hdr.h".)
  We use Grid macros to declare and implement specific instances of Grid
*/

implement_Grid2(double)

/*
  For each grid type we must have a declaration
  and corresponding implementation.
  Here we implement a 2-D Grid class of double
  The corresponding declaration appears before it in the header file "hdr.h"
  Each module using a particular Grid type must declare it, i.e. in "hdr.h"
  The implementation must appear in exactly one .C module (generally main)
  This is a standard C++ way of emulating templates. 
  Implementations look like function calls, but in fact they invoke macros
  Because type double is a built-in type, we do no have to 
  define any pack-unpack routines for double.
  However, if the type were a user-defined class or structure,
  then we would need to provide hints to LPARX about how to
  gather and scatter the data type. This is relatively easy
  to do and is described in the User's Guide.
*/

//
//  forward and external procedure declarations
//

void initData(Grid2(double)& U);
double relax(Grid2(double)& U, const double& h);

//
// Beginning of main() -- define defaults for the problem
//

#define M 40
#define N 50
#define MAX_ITERATIONS 10
#define DELTA_RESIDUAL 1.0e-2

int main(const int argc, char **argv)
  {
   /*
     Define the interior problem domain to be [1:M, 1:N]
     A Region2 is a first-class object that describes an index domain.
   */

   Region2 Interior(1, 1, M, N);
   Region2 W = grow(Interior, 1);
   Grid2(double) U(W);
   const double h = 1.0;

   // Set up the initial conditions

   initData(U);

   // Perform the iterations

   for (int i = 0; i < MAX_ITERATIONS; i++)
     {
      const double dU = relax(U, h);
      cout << "Delta at iteration " << i << ": " << dU << endl;
      if (dU < DELTA_RESIDUAL) break;
     }
  }