relax.c

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

#include "redblack2d.h"

//
// Fortran name mangling is architecture-specific, so define a mangled
// name appropriate to your architecture
//

#if (defined(nCUBE2) || defined(CRAYC90))
#define redblack5 REDBLACK5
#define residual5 RESIDUAL5
#else
#define redblack5 redblack5_
#define residual5 residual5_
#endif

extern "C" {
void redblack5(double *u,
               const int *const nl, const int *const nh,
               const double *h, const int *const offset, double *du);

}

extern void fillpatch(XArray1(Grid2(double))& U);

double relax(XArray1(Grid2(double))& U, const double& h)
  {
   /*	
      Call fillpatch() to acquire ghost cell values.
   */

   fillpatch(U);

   /*
     Call a Fortran routine to perform relaxation on each grid.

     The for_all_1 parallel loop is described in main(). 

     This is an example of the C++ - Fortran interface. The C++
     grid class uses Fortran-style array ordering of the indices.
     So, to call Fortran, we need to pass the beginning of the
     grid array and also the x and y dimensions of the array. 

     The member function X.extents(Y), where X is a Region,
     returns the size of the Yth dimension of the Region object
     X. The member function X.data(), where X is a grid,
     returns a pointer the buffer holding the grid data of object X.

     Remember the parameter passing conventions in Fortran --
     use pointers to scalars. redblack5 is a Fortran subroutine
     that performs a red-black jacobi update on a 2D array
     of REAL*8 values (see file redblack5.F).
   */

   const int red = 0, black = 1;
   double r0=0.0, r1=0.0;
   for_all_1(i, U)
      Region2 inside = grow(U(i).region(), -1);
      const Point2 nl = inside.lower();
      const Point2 nh = inside.upper();
      redblack5(U(i).data(), nl, nh, &h, &red, &r0);
      redblack5(U(i).data(), nl, nh, &h, &black, &r1);
      if (r1 > r0) r0 = r1;
   end_for_all_no_sync

   mpReduceMax(&r0);
   return(r0);
  }