Yet more stuff. Nw have routines to create pGAs.

Tom


git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGH/trunk@131 b61c5cb5-eaca-4651-9a7a-d64986f99364

1 files changed, 653 insertions, 354 deletions

diff --git a/src/SetupPGV.c b/src/SetupPGV.c
index 8b1c676..c41ea1e 100644
--- a/src/SetupPGV.c
+++ b/src/SetupPGV.c
@@ -12,17 +12,21 @@ static char *rcsid = "$Header$";
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include <string.h>
 
 /* Stuff to do standalone test */
 #define TEST_SETUPPGV
 
 #ifdef TEST_SETUPPGV
-#define PUGH_NSTAGGER 4
+#define PUGH_NSTAGGER   4
 #define PUGH_NO_STAGGER 0
-#define PUGH_STAGGER 1
-#define PUGH_VERTEXCTR 0
+#define PUGH_STAGGER    1
+#define PUGH_VERTEXCTR  0
 #define PUGH_FACECTR(i) i+1
-#define CCTK_REAL double
+#define PUGH_NOSTORAGE  0
+#define PUGH_STORAGE    1
+#define PUGH_NOCOMM     0
+#define PUGH_COMM       1
 #endif /* TEST_SETUPPGV */
 
 typedef enum {pgv_none, pgv_scalar, pgv_array, pgv_gf} pgv_type;
@@ -43,68 +47,69 @@ typedef struct PGS
 
 typedef struct PGExtras
 {
-  int   dim;                    /* dimension of GF */
+  int	dim;			/* dimension of GA */
 
-  int *nsize;                  /* The global size of the array */
+  int *nsize;		       /* The global size of the array */
 
   /* Processor group layouts */
-  CCTK_REAL maxskew;		/* Maximum point skew */
-  int     **lb;			/* Lower bound (nprocs X dim) for each proc */
-  int     **ub;			/* Upper bound (same sizes) */
-  int      *lnsize;		/* Size on this processor */
-  int       npoints;		/* LOCAL number of points on this proc. */
-  int      *rnpoints;		/* Number of points on each proc */
-  int     **rnsize;	        /* [#points on a proc][in each dir] */
+  double   maxskew;		/* Maximum point skew */
+  int	  **lb;			/* Lower bound (nprocs X dim) for each proc */
+  int	  **ub;			/* Upper bound (same sizes) */
+  int	   *lnsize;		/* Size on this processor */
+  int	    npoints;		/* LOCAL number of points on this proc. */
+  int	   *rnpoints;		/* Number of points on each proc */
+  int	  **rnsize;		/* [#points on a proc][in each dir] */
 
   /* Ghosts and overlaps. */
-  int      *nghostzones;	/* Width of ghost zone */
-  int      *ownership[PUGH_NSTAGGER][2];	
-                                /* The box owned in each direction. */
-                                /* [stagger][min/max][dir] */
+  int	   *nghostzones;	/* Width of ghost zone */
+  int	   *ownership[PUGH_NSTAGGER][2];	
+				/* The box owned in each direction. */
+				/* [stagger][min/max][dir] */
 				
-  int     **ghosts[PUGH_NSTAGGER][2];	
-                                /* The ghost zones on each face. */
-                                /* [stagger][min/max][face][dir] */
+  int	  **ghosts[PUGH_NSTAGGER][2];	
+				/* The ghost zones on each face. */
+				/* [stagger][min/max][face][dir] */
 				
-  int     **overlap[PUGH_NSTAGGER][2];	
-                                /* The overlap region owned on each face. */
-			        /* [stagger][min/max][face][ijk] */
+  int	  **overlap[PUGH_NSTAGGER][2];	
+				/* The overlap region owned on each face. */
+				/* [stagger][min/max][face][ijk] */
 
 } pGExtras;
 
 typedef struct PGA
 {
-  char  *name;               /* The name of the grid function */
-  int   gfno;                /* My ID number in my GH parent. */
-  void  *padddata;           /* Storage for the data. */
-  void  *data;               /* See the note above. */
-  /*     buffer_sz[2*dim] */
-  int   *buffer_sz;        /* Size of the face ghost zones */
-  void  **send_buffer;       /* Storage for buffered Comm if */
-  void  **recv_buffer;       /* we don't use derived types */
-  int   commflag;            /* What is the comm flag set to? */
-  /*     do_comm[2*dim] */
-  int   *docomm;           /* Do we do comm or not? */
-  int   storage;             /* Do we have storage or not? */
-  int   stagger;             /* Only Vertex Centered now... */
-
-  struct PGH      *parentGH; /* The GH to which I belong */
-                             /* Note this is struct PGH whic is
-                                typedeffed to pGH in pGH.h, but
-                                that is included AFTER this so we 
-                                need the full name for the lookahead
-                                structure reference thingy.
-                             */
-
-  int   varsize;             /* The size of the data */
-  int   vtype;               /* The type of the data */
+  char	*name;		     /* The name of the grid function */
+  int	id;		   /* My ID number in my GH parent. */
+  void	*padddata;	     /* Storage for the data. */
+  void	*data;		     /* See the note above. */
+  /*	 buffer_sz[2*dim] */
+  int	*buffer_sz;	   /* Size of the face ghost zones */
+  void	**send_buffer;	     /* Storage for buffered Comm if */
+  void	**recv_buffer;	     /* we don't use derived types */
+  int	commflag;	     /* What is the comm flag set to? */
+  /*	 do_comm[2*dim] */
+  int	*docomm;	   /* Do we do comm or not? */
+  int	storage;	     /* Do we have storage or not? */
+  int	stagger;	     /* Only Vertex Centered now... */
+
+  void	*parent;	     /* The GH to which I belong */
+			     /* Note this is struct PGH whic is
+				typedeffed to pGH in pGH.h, but
+				that is included AFTER this so we 
+				need the full name for the lookahead
+				structure reference thingy.
+			     */
+
+  int	varsize;	     /* The size of the data */
+  int	vtype;		     /* The type of the data */
 
 #ifdef MPI
-  MPI_Request sreq[6], rreq[6];  /* Comm requests and statuses. */
+  MPI_Request *sreq, *rreq;	 /* Comm requests and statuses. */
   MPI_Status ms;
 #endif
 
   pGExtras *extras;
+  pConnectivity *connectivity;
 } pGA;
 
 typedef struct PGV
@@ -116,13 +121,13 @@ typedef struct PGV
 } pGV;
 
 pGExtras *pugh_SetupPGExtras(int dim, 
-                             int periodic,
-                             int staggertype,
-                             int *sh,
-                             int *nghosts,
-                             int total_procs,
-                             int *nprocs,
-                             int this_proc)
+			     int periodic,
+			     int stagger,
+			     int *sh,
+			     int *nghosts,
+			     int total_procs,
+			     int *nprocs,
+			     int this_proc)
 {
   int error;
   pGExtras *this;
@@ -133,17 +138,17 @@ pGExtras *pugh_SetupPGExtras(int dim,
   if(this)
   {
     error = pugh_SetupPGExtrasMemory(dim, 
-                                     total_procs,
-                                     nprocs,
-                                     this);
+				     total_procs,
+				     nprocs,
+				     this);
 
     if(!error)
     {
-      pugh_SetupPGExtrasSizes   (dim, periodic, staggertype, sh, nghosts, 
-                                 total_procs, nprocs, this_proc,this);
-      pugh_SetupPGExtrasOwnership(dim, periodic, staggertype, sh, nghosts, 
-                                 total_procs, nprocs, this_proc, this);
-    }      
+      pugh_SetupPGExtrasSizes	(dim, periodic, stagger, sh, nghosts, 
+				 total_procs, nprocs, this_proc,this);
+      pugh_SetupPGExtrasOwnership(dim, periodic, stagger, sh, nghosts, 
+				 total_procs, nprocs, this_proc, this);
+    }	   
 
   }
 
@@ -152,7 +157,7 @@ pGExtras *pugh_SetupPGExtras(int dim,
 
  /*@@
    @routine    pugh_SetupConnectivity
-   @date       Fri Nov  5 11:32:12 1999
+   @date       Fri Nov	5 11:32:12 1999
    @author     Tom Goodale
    @desc 
    Create a connectivity struvcture containing 
@@ -167,9 +172,9 @@ pGExtras *pugh_SetupPGExtras(int dim,
 
 @@*/
 pConnectivity *pugh_SetupConnectivity(int dim, 
-                                      int total_procs, 
-                                      int *nprocs,
-                                      int periodic)
+				      int total_procs, 
+				      int *nprocs,
+				      int periodic)
 {
   pConnectivity *this;
 
@@ -189,19 +194,19 @@ pConnectivity *pugh_SetupConnectivity(int dim,
 
       if(this->neighbours[0])
       {
-        for(i = 1; i < total_procs; i++)
-        {
-          this->neighbours[i] = this->neighbours[0]+(2*dim*i);
-        }
+	for(i = 1; i < total_procs; i++)
+	{
+	  this->neighbours[i] = this->neighbours[0]+(2*dim*i);
+	}
 
-        this->periodic = periodic;
+	this->periodic = periodic;
       }
       else
       {
-        free(this->nprocs);
-        free(this->neighbours);
-        free(this);
-        this = NULL;
+	free(this->nprocs);
+	free(this->neighbours);
+	free(this);
+	this = NULL;
       }
     }
     else
@@ -232,7 +237,7 @@ pConnectivity *pugh_SetupConnectivity(int dim,
 
  /*@@
    @routine    pugh_GenerateTopology
-   @date       Fri Nov  5 11:31:21 1999
+   @date       Fri Nov	5 11:31:21 1999
    @author     Tom Goodale
    @desc 
    Generate the appropriate processor topology for this processor
@@ -266,11 +271,11 @@ int pugh_GenerateTopology(int dim, int total_procs, int *nprocs)
       free_dims--;
       if(used_procs) 
       {
-        used_procs *= nprocs[i];
+	used_procs *= nprocs[i];
       }
       else
       {
-        used_procs = nprocs[i];
+	used_procs = nprocs[i];
       }
     }
   }
@@ -301,12 +306,12 @@ int pugh_GenerateTopology(int dim, int total_procs, int *nprocs)
     {
       if(!nprocs[i])
       {
-        nprocs[i] = pow(free_procs, 1.0/free_dims);
+	nprocs[i] = pow(free_procs, 1.0/free_dims);
 
-        while(free_procs % nprocs[i]) nprocs[i]--;
+	while(free_procs % nprocs[i]) nprocs[i]--;
 
-        free_procs /= nprocs[i];
-        free_dims--;
+	free_procs /= nprocs[i];
+	free_dims--;
       }
     }
   }
@@ -316,7 +321,7 @@ int pugh_GenerateTopology(int dim, int total_procs, int *nprocs)
 
  /*@@
    @routine    pugh_GenerateNeighbours
-   @date       Mon Nov  8 08:15:08 1999
+   @date       Mon Nov	8 08:15:08 1999
    @author     Tom Goodale
    @desc 
    Works out the array of neighbouring processors for
@@ -333,6 +338,7 @@ int pugh_GenerateNeighbours(int dim, int total_procs, int *nprocs, int **neighbo
 {
   int retval;
   int i;
+
   int idim;
   int *pos;
   int temp;
@@ -347,39 +353,39 @@ int pugh_GenerateNeighbours(int dim, int total_procs, int *nprocs, int **neighbo
       
       for(idim = 0; idim < dim ; idim++)
       {
-        /* Deal with minus neighbour in this direction */
-        pos[idim]--;
-        
-        if(pos[idim] > -1)
-        {
-          neighbours[i][idim*2] = pugh_ComposeIJK(dim, nprocs, pos);
-        }
-        else if(periodic)
-        {
-          temp = pos[idim];
-          pos[idim] = nprocs[idim]-1;
-          neighbours[i][idim*2] = pugh_ComposeIJK(dim, nprocs, pos);
-          pos[idim] = temp;
-        }
+	/* Deal with minus neighbour in this direction */
+	pos[idim]--;
+	
+	if(pos[idim] > -1)
+	{
+	  neighbours[i][idim*2] = pugh_ComposeIJK(dim, nprocs, pos);
+	}
+	else if(periodic)
+	{
+	  temp = pos[idim];
+	  pos[idim] = nprocs[idim]-1;
+	  neighbours[i][idim*2] = pugh_ComposeIJK(dim, nprocs, pos);
+	  pos[idim] = temp;
+	}
 
-        pos[idim]++;
-        
-        /* Deal with plus neighbour in this direction */
-        pos[idim]++;
+	pos[idim]++;
+	
+	/* Deal with plus neighbour in this direction */
+	pos[idim]++;
 
-        if(pos[idim] < nprocs[idim])
-        {
-          neighbours[i][idim*2+1] = pugh_ComposeIJK(dim, nprocs, pos);
-        }
-        else if(periodic)
-        {
-          temp = pos[idim];
-          pos[idim] = 0;
-          neighbours[i][idim*2+1] = pugh_ComposeIJK(dim, nprocs, pos);
-          pos[idim] = temp;
-        }
-        
-        pos[idim]--;
+	if(pos[idim] < nprocs[idim])
+	{
+	  neighbours[i][idim*2+1] = pugh_ComposeIJK(dim, nprocs, pos);
+	}
+	else if(periodic)
+	{
+	  temp = pos[idim];
+	  pos[idim] = 0;
+	  neighbours[i][idim*2+1] = pugh_ComposeIJK(dim, nprocs, pos);
+	  pos[idim] = temp;
+	}
+	
+	pos[idim]--;
       }
     }
 
@@ -399,7 +405,7 @@ int pugh_GenerateNeighbours(int dim, int total_procs, int *nprocs, int **neighbo
 
  /*@@
    @routine    pugh_DecomposeIJK
-   @date       Fri Nov  5 11:29:43 1999
+   @date       Fri Nov	5 11:29:43 1999
    @author     Tom Goodale
    @desc 
    Decompose an ijk index into seperate components.
@@ -420,7 +426,7 @@ int pugh_DecomposeIJK(int dim, int ijk, int *nprocs, int *pos)
    *
    * e.g. ijk = i+nx*(j+ny*k)
    * => i = ijk % nx
-   *    ijk/nx = j + ny*k
+   *	ijk/nx = j + ny*k
    */
 
   for(idim = 0; idim < dim; idim++)
@@ -435,7 +441,7 @@ int pugh_DecomposeIJK(int dim, int ijk, int *nprocs, int *pos)
 
  /*@@
    @routine    pugh_ComposeIJK
-   @date       Fri Nov  5 11:29:43 1999
+   @date       Fri Nov	5 11:29:43 1999
    @author     Tom Goodale
    @desc 
    Compose an ijk index from seperate components.
@@ -449,8 +455,8 @@ int pugh_DecomposeIJK(int dim, int ijk, int *nprocs, int *pos)
 
 @@*/
 int pugh_ComposeIJK(int dim, 
-                    int *nprocs, 
-                    int *pos)
+		    int *nprocs, 
+		    int *pos)
 {
   int ijk;
   int idim;
@@ -472,7 +478,7 @@ int pugh_ComposeIJK(int dim,
 
  /*@@
    @routine    pugh_SetupPGExtrasMemory
-   @date       Mon Nov  8 08:16:02 1999
+   @date       Mon Nov	8 08:16:02 1999
    @author     Tom Goodale
    @desc 
    Allocate memory for the members of the pGExtras structure.
@@ -485,9 +491,9 @@ int pugh_ComposeIJK(int dim,
 
 @@*/
 int pugh_SetupPGExtrasMemory(int dim, 
-                             int total_procs,
-                             int *nprocs,
-                             pGExtras *this)
+			     int total_procs,
+			     int *nprocs,
+			     pGExtras *this)
 {
   int retcode;
   int i,j,k;
@@ -499,67 +505,67 @@ int pugh_SetupPGExtrasMemory(int dim,
     /* Do it in stages.
      * First: things depending on the number of processors 
      */
-    this->lb       = (int **)malloc(total_procs*sizeof(int *));
-    this->ub       = (int **)malloc(total_procs*sizeof(int *));
+    this->lb	   = (int **)malloc(total_procs*sizeof(int *));
+    this->ub	   = (int **)malloc(total_procs*sizeof(int *));
     this->rnsize   = (int **)malloc(total_procs*sizeof(int *));
     this->rnpoints = (int *) malloc(total_procs*sizeof(int));
 
     /* Things just depending on dimension */
     this->nghostzones = (int *)malloc(dim*sizeof(int));
-    this->nsize       = (int *)malloc(dim*sizeof(int));
+    this->nsize	      = (int *)malloc(dim*sizeof(int));
     this->lnsize      = (int *)malloc(dim*sizeof(int));
 
     /* Check all the above succeeded and then get memory for 
      * arrays hanging off the above. 
      */
-    if(this->lb          &&
-       this->ub          &&
-       this->rnsize      &&
-       this->rnpoints    &&
+    if(this->lb		 &&
+       this->ub		 &&
+       this->rnsize	 &&
+       this->rnpoints	 &&
        this->nghostzones &&
-       this->nsize       &&
+       this->nsize	 &&
        this->lnsize)
     {
       this->lb[0]     = (int *)malloc(total_procs  *dim*sizeof(int));
       this->ub[0]     = (int *)malloc(total_procs  *dim*sizeof(int));
       this->rnsize[0] = (int *)malloc(total_procs*2*dim*sizeof(int));
 
-      if(this->lb[0]     &&
-         this->ub[0]     &&
-         this->rnsize[0])
+      if(this->lb[0]	 &&
+	 this->ub[0]	 &&
+	 this->rnsize[0])
       {
-        for (i = 1; i < total_procs; i++)
-        {
-          this->lb[i]     = this->lb[0]     + i*dim;
-          this->ub[i]     = this->ub[0]     + i*dim;
-          this->rnsize[i] = this->rnsize[0] + i*2*dim;
-        }
+	for (i = 1; i < total_procs; i++)
+	{
+	  this->lb[i]	  = this->lb[0]	    + i*dim;
+	  this->ub[i]	  = this->ub[0]	    + i*dim;
+	  this->rnsize[i] = this->rnsize[0] + i*2*dim;
+	}
       }
       else
       {
-        /* Free inner arrays */
-        free(this->lb[0]);
-        this->lb[0] = NULL;
-        free(this->ub[0]);
-        this->ub[0] = NULL;
-        free(this->rnsize[0]);
-        this->rnsize[0] = NULL;
-
-        /* Free toplevel arrays */
-        free(this->lb);
-        this->lb = NULL;
-        free(this->ub);
-        this->ub = NULL;
-        free(this->rnsize);
-        this->rnsize = NULL;
-        free(this->rnpoints);
-        this->rnpoints = NULL;
-        free(this->nghostzones);
-        this->nghostzones = NULL;               
-        free(this->nsize);
-        this->nsize = NULL;               
-        free(this->lnsize);
-        this->lnsize = NULL;               
+	/* Free inner arrays */
+	free(this->lb[0]);
+	this->lb[0] = NULL;
+	free(this->ub[0]);
+	this->ub[0] = NULL;
+	free(this->rnsize[0]);
+	this->rnsize[0] = NULL;
+
+	/* Free toplevel arrays */
+	free(this->lb);
+	this->lb = NULL;
+	free(this->ub);
+	this->ub = NULL;
+	free(this->rnsize);
+	this->rnsize = NULL;
+	free(this->rnpoints);
+	this->rnpoints = NULL;
+	free(this->nghostzones);
+	this->nghostzones = NULL;		
+	free(this->nsize);
+	this->nsize = NULL;		  
+	free(this->lnsize);
+	this->lnsize = NULL;		   
       }
     }
     else
@@ -574,122 +580,122 @@ int pugh_SetupPGExtrasMemory(int dim,
       free(this->rnpoints);
       this->rnpoints = NULL;
       free(this->nghostzones);
-      this->nghostzones = NULL;               
+      this->nghostzones = NULL;		      
       free(this->nsize);
-      this->nsize = NULL;               
+      this->nsize = NULL;		
       free(this->lnsize);
-      this->lnsize = NULL;               
+      this->lnsize = NULL;		 
     }
       
 
-    if(this->lb          &&
-       this->ub          &&
-       this->rnsize      &&
-       this->rnpoints    &&
+    if(this->lb		 &&
+       this->ub		 &&
+       this->rnsize	 &&
+       this->rnpoints	 &&
        this->nghostzones &&
-       this->nsize       &&
+       this->nsize	 &&
        this->lnsize)
     {
       retcode = 0;
       for (i = 0 ; i < PUGH_NSTAGGER; i++)
       {
-        for (j = 0; j < 2; j++)
-        {
-          this->ownership[i][j]    = (int *) malloc(dim*sizeof(int));
-          this->ghosts[i][j]       = (int **)malloc(2*dim*sizeof(int *));
-          this->overlap[i][j]      = (int **)malloc(2*dim*sizeof(int *));
-          if(this->ghosts[i][j] &&
-             this->overlap[i][j])
-          {
-            this->ghosts[i][j][0]  = (int *)malloc(2*dim*dim*sizeof(int));
-            this->overlap[i][j][0] = (int *)malloc(2*dim*dim*sizeof(int));
-            for (k=1; k < 2*dim; k++)
-            {
-              this->ghosts[i][j][k]  = this->ghosts[i][j][0] + k*dim;
-              this->overlap[i][j][k] = this->overlap[i][j][0] + k*dim;
-            }
-          }
-          else
-          {
-            free(this->ownership[i][j]);
-            this->ownership[i][j] = NULL;
-
-            free(this->ghosts[i][j]);
-            this->ghosts[i][j] = NULL;
-
-            free(this->overlap[i][j]);
-            this->overlap[i][j] = NULL;
-            retcode = 1;
-            break;
-          }          
-        }
-        if(retcode)
-        {
-          for(j=1; j >=0 ; j--)
-          {
-              free(this->ownership[i][j]);
-              this->ownership[i][j] = NULL;
-              
-              free(this->ghosts[i][j]);
-              this->ghosts[i][j] = NULL;
-              
-              free(this->overlap[i][j]);
-              this->overlap[i][j] = NULL;
-              
-          }
-          break;
-        }
+	for (j = 0; j < 2; j++)
+	{
+	  this->ownership[i][j]	   = (int *) malloc(dim*sizeof(int));
+	  this->ghosts[i][j]	   = (int **)malloc(2*dim*sizeof(int *));
+	  this->overlap[i][j]	   = (int **)malloc(2*dim*sizeof(int *));
+	  if(this->ghosts[i][j] &&
+	     this->overlap[i][j])
+	  {
+	    this->ghosts[i][j][0]  = (int *)malloc(2*dim*dim*sizeof(int));
+	    this->overlap[i][j][0] = (int *)malloc(2*dim*dim*sizeof(int));
+	    for (k=1; k < 2*dim; k++)
+	    {
+	      this->ghosts[i][j][k]  = this->ghosts[i][j][0] + k*dim;
+	      this->overlap[i][j][k] = this->overlap[i][j][0] + k*dim;
+	    }
+	  }
+	  else
+	  {
+	    free(this->ownership[i][j]);
+	    this->ownership[i][j] = NULL;
+
+	    free(this->ghosts[i][j]);
+	    this->ghosts[i][j] = NULL;
+
+	    free(this->overlap[i][j]);
+	    this->overlap[i][j] = NULL;
+	    retcode = 1;
+	    break;
+	  }	     
+	}
+	if(retcode)
+	{
+	  for(j=1; j >=0 ; j--)
+	  {
+	      free(this->ownership[i][j]);
+	      this->ownership[i][j] = NULL;
+	      
+	      free(this->ghosts[i][j]);
+	      this->ghosts[i][j] = NULL;
+	      
+	      free(this->overlap[i][j]);
+	      this->overlap[i][j] = NULL;
+	      
+	  }
+	  break;
+	}
       }
       if(retcode)
       {
-        /* Loop back through the arrays freeing things */
-        for(i--; i >=0; i--)
-        {
-          for(j=1; j >=0 ; j--)
-          {
-            free(this->ghosts[i][j][0]);
-            free(this->overlap[i][j][0]);
-
-            free(this->ownership[i][j]);
-            this->ownership[i][j] = NULL;
-            
-            free(this->ghosts[i][j]);
-            this->ghosts[i][j] = NULL;
-          
-            free(this->overlap[i][j]);
-            this->overlap[i][j] = NULL;
-          }
+	/* Loop back through the arrays freeing things */
+	for(i--; i >=0; i--)
+	{
+	  for(j=1; j >=0 ; j--)
+	  {
+	    free(this->ghosts[i][j][0]);
+	    free(this->overlap[i][j][0]);
+
+	    free(this->ownership[i][j]);
+	    this->ownership[i][j] = NULL;
+	    
+	    free(this->ghosts[i][j]);
+	    this->ghosts[i][j] = NULL;
+	  
+	    free(this->overlap[i][j]);
+	    this->overlap[i][j] = NULL;
+	  }
 
-          free(this->ownership[i]);
-          free(this->ghosts[i]);
-          free(this->overlap[i]);
-        }
+	  free(this->ownership[i]);
+	  free(this->ghosts[i]);
+	  free(this->overlap[i]);
+	}
 
-        /* Free the stuff originally allocated */
-
-        free(this->lb[0]);
-        this->lb[0] = NULL;
-        free(this->ub[0]);
-        this->ub[0] = NULL;
-        free(this->rnsize[0]);
-        this->rnsize[0] = NULL;
-
-        free(this->lb);
-        this->lb = NULL;
-        free(this->ub);
-        this->ub = NULL;
-        free(this->rnsize);
-        this->rnsize = NULL;
-        free(this->rnpoints);
-        this->rnpoints = NULL;
-        free(this->nghostzones);
-        this->nghostzones = NULL;               
-        free(this->nsize);
-        this->nsize = NULL;               
-        free(this->lnsize);
-        this->lnsize = NULL;               
+	/* Free the stuff originally allocated */
+
+	free(this->lb[0]);
+	this->lb[0] = NULL;
+	free(this->ub[0]);
+	this->ub[0] = NULL;
+	free(this->rnsize[0]);
+	this->rnsize[0] = NULL;
+
+	free(this->lb);
+	this->lb = NULL;
+	free(this->ub);
+	this->ub = NULL;
+	free(this->rnsize);
+	this->rnsize = NULL;
+	free(this->rnpoints);
+	this->rnpoints = NULL;
+	free(this->nghostzones);
+	this->nghostzones = NULL;		
+	free(this->nsize);
+	this->nsize = NULL;		  
+	free(this->lnsize);
+	this->lnsize = NULL;		   
       }
-    }          
+    }	       
   }
   else
   {
@@ -701,7 +707,7 @@ int pugh_SetupPGExtrasMemory(int dim,
 
  /*@@
    @routine    pugh_SetupPGExtrasSizes
-   @date       Mon Nov  8 08:59:33 1999
+   @date       Mon Nov	8 08:59:33 1999
    @author     Tom Goodale
    @desc 
    Sets up the size information in the pGExtras 
@@ -715,14 +721,14 @@ int pugh_SetupPGExtrasMemory(int dim,
 
 @@*/
 int pugh_SetupPGExtrasSizes(int dim, 
-                            int periodic,
-                            int staggertype,
-                            int *sh,
-                            int *nghosts,
-                            int total_procs,
-                            int *nprocs,
-                            int this_proc,
-                            pGExtras *this)
+			    int periodic,
+			    int stagger,
+			    int *sh,
+			    int *nghosts,
+			    int total_procs,
+			    int *nprocs,
+			    int this_proc,
+			    pGExtras *this)
 {
   int dir;
 
@@ -736,7 +742,7 @@ int pugh_SetupPGExtrasSizes(int dim,
 	(-sh[dir] - 2*nghosts[dir]) +
 	2 * (-sh[dir] - nghosts[dir]);
 
-      if (staggertype == PUGH_STAGGER) 
+      if (stagger == PUGH_STAGGER) 
       {
 	this->nsize[dir] -= nprocs[dir]-1;
       }
@@ -755,11 +761,11 @@ int pugh_SetupPGExtrasSizes(int dim,
 
   /* Setup the bounding box stuff */
 
-  pugh_SetupBoundingBox(dim, periodic, staggertype, sh, nghosts, total_procs, nprocs, this); 
+  pugh_SetupBoundingBox(dim, periodic, stagger, sh, nghosts, total_procs, nprocs, this); 
 
   /* Set the remote sizes */
 
-  pugh_SetupRemoteSizes(dim, periodic, staggertype, sh, nghosts, total_procs, nprocs, this); 
+  pugh_SetupRemoteSizes(dim, periodic, stagger, sh, nghosts, total_procs, nprocs, this); 
   
   /* Set the local sizes */
 
@@ -768,6 +774,7 @@ int pugh_SetupPGExtrasSizes(int dim,
     this->lnsize[dir] = this->rnsize[this_proc][dir];
   }
 
+  this->npoints = this->rnpoints[this_proc];
 
   return 0;
 }
@@ -775,7 +782,7 @@ int pugh_SetupPGExtrasSizes(int dim,
 
  /*@@
    @routine    pugh_SetupPGExtrasOwnership
-   @date       Mon Nov  8 09:00:10 1999
+   @date       Mon Nov	8 09:00:10 1999
    @author     Tom Goodale
    @desc 
    Sets up ownership, overlap, ghostzones, etc 
@@ -791,14 +798,14 @@ int pugh_SetupPGExtrasSizes(int dim,
 
 @@*/
 int pugh_SetupPGExtrasOwnership(int dim, 
-                               int periodic,
-                               int staggertype,
-                               int *sh,
-                               int *nghosts,
-                               int total_procs,
-                               int *nprocs,
-                               int this_proc,
-                               pGExtras *this)
+			       int periodic,
+			       int stagger,
+			       int *sh,
+			       int *nghosts,
+			       int total_procs,
+			       int *nprocs,
+			       int this_proc,
+			       pGExtras *this)
 {
   int tmp;
   int dir, idir;
@@ -814,10 +821,10 @@ int pugh_SetupPGExtrasOwnership(int dim,
   for (dir = 0 ; dir < dim; dir++)
   {
     this->ownership[PUGH_VERTEXCTR][0][dir] = (this->lb[this_proc][dir] == 0 ? 
-                                               0 : this->nghostzones[dir]);
+					       0 : this->nghostzones[dir]);
     this->ownership[PUGH_VERTEXCTR][1][dir]=(this->ub[this_proc][dir] == this->nsize[dir]-1 ? 
-                                             this->lnsize[dir] : this->lnsize[dir] - 
-                                             this->nghostzones[dir]);
+					     this->lnsize[dir] : this->lnsize[dir] - 
+					     this->nghostzones[dir]);
   }
 
   /* correct for periodic identification : Tue Jun 17 08:40:15 CDT 1997 */
@@ -866,7 +873,7 @@ int pugh_SetupPGExtrasOwnership(int dim,
 	   Luckily iend is aready set to nx ... */
 	tmp = iend;
 	istart = tmp - 2 * this->nghostzones[idir];
-	iend   = tmp -     this->nghostzones[idir];
+	iend   = tmp -	   this->nghostzones[idir];
       }
 
       this->overlap[PUGH_VERTEXCTR][0][dir][idir] = istart;
@@ -896,8 +903,8 @@ int pugh_SetupPGExtrasOwnership(int dim,
       if (dir == 2*idir+1) 
       {
       /* Want to go from nx - 2 sw to nx - sw
-         (eg nx-2 to < nx-1, nx-4 to < nx-2)
-         Luckily iend is aready set to nx ... */
+	 (eg nx-2 to < nx-1, nx-4 to < nx-2)
+	 Luckily iend is aready set to nx ... */
 	tmp = iend;
 	istart = tmp - this->nghostzones[idir];
 	iend   = tmp;
@@ -911,7 +918,7 @@ int pugh_SetupPGExtrasOwnership(int dim,
   }
 
   /* if staggering, set up ownership, overlaps, ghosts */
-  if(staggertype == PUGH_STAGGER) 
+  if(stagger == PUGH_STAGGER) 
   {
     /* copy ownership from PUGH_VERTEXCTR */
     for(i = 1; i < PUGH_NSTAGGER; i++)
@@ -919,9 +926,9 @@ int pugh_SetupPGExtrasOwnership(int dim,
       for(j = 0; j < dim; j++)
       {
 	for(k = 0; k < 2; k++)
-        {
+	{
 	  this->ownership[i][k][j] = this->ownership[PUGH_VERTEXCTR][k][j];
-        }
+	}
       }
     }
 
@@ -931,16 +938,16 @@ int pugh_SetupPGExtrasOwnership(int dim,
       for (idir = 0; idir < dim; idir++)
       {	
 	if(this->ub[this_proc][idir] != this->nsize[idir]-1)
-        {
+	{
 	  this->ownership[i][1][idir] --;
-        }
+	}
 	else
-        {
+	{
 	  if(i == PUGH_FACECTR(idir))
-          {
+	  {
 	    this->ownership[i][1][idir] --;
-          }
-        }
+	  }
+	}
       }
     }
 
@@ -949,33 +956,33 @@ int pugh_SetupPGExtrasOwnership(int dim,
     {
       for(dir = 0; dir < 2*dim; dir++)
       {
-        for (idir = 0; idir < dim; idir++)
-        {
-          istart = 0;  
-          iend   = this->lnsize[idir];
-          if(staglcv == PUGH_FACECTR(idir))
-          {
-            iend --;
-          }
-          if(dir == 2*idir)
-          {
-            istart =     this->nghostzones[idir];
-            iend   = 2 * this->nghostzones[idir];
-            if(staglcv != PUGH_FACECTR(idir))
-            {
-              istart ++;
-              iend ++;
-            }
-          }
-          if(dir == 2*idir+1)
-          {
-            istart = this->lnsize[idir] - 2 * this->nghostzones[idir] -1;
-            iend   = istart                 + this->nghostzones[idir];
-          }
-          
-          this->overlap[staglcv][0][dir][idir] = istart;
-          this->overlap[staglcv][1][dir][idir] = iend;
-        }
+	for (idir = 0; idir < dim; idir++)
+	{
+	  istart = 0;  
+	  iend	 = this->lnsize[idir];
+	  if(staglcv == PUGH_FACECTR(idir))
+	  {
+	    iend --;
+	  }
+	  if(dir == 2*idir)
+	  {
+	    istart =	 this->nghostzones[idir];
+	    iend   = 2 * this->nghostzones[idir];
+	    if(staglcv != PUGH_FACECTR(idir))
+	    {
+	      istart ++;
+	      iend ++;
+	    }
+	  }
+	  if(dir == 2*idir+1)
+	  {
+	    istart = this->lnsize[idir] - 2 * this->nghostzones[idir] -1;
+	    iend   = istart		    + this->nghostzones[idir];
+	  }
+	  
+	  this->overlap[staglcv][0][dir][idir] = istart;
+	  this->overlap[staglcv][1][dir][idir] = iend;
+	}
       }
     }
        
@@ -987,11 +994,11 @@ int pugh_SetupPGExtrasOwnership(int dim,
 	for(idir=0;idir<dim;idir++)
 	{
 	  istart = 0;  
-	  iend   = this->lnsize[idir];
+	  iend	 = this->lnsize[idir];
 	  if(staglcv == PUGH_FACECTR(idir))
-          {
-            iend --;
-          }
+	  {
+	    iend --;
+	  }
 	  if(dir == 2*idir)
 	  {
 	    istart = 0;
@@ -1019,7 +1026,7 @@ int pugh_SetupPGExtrasOwnership(int dim,
 
  /*@@
    @routine    pugh_SetupBoundingBox
-   @date       Mon Nov  8 09:03:40 1999
+   @date       Mon Nov	8 09:03:40 1999
    @author     Tom Goodale
    @desc 
    Sets up the bounding box info for a pgExtras structure.
@@ -1032,13 +1039,13 @@ int pugh_SetupPGExtrasOwnership(int dim,
 
 @@*/
 int pugh_SetupBoundingBox(int dim, 
-                          int periodic,
-                          int staggertype,
-                          int *sh,
-                          int *nghosts,
-                          int total_procs,
-                          int *nprocs,
-                          pGExtras *this)
+			  int periodic,
+			  int stagger,
+			  int *sh,
+			  int *nghosts,
+			  int total_procs,
+			  int *nprocs,
+			  pGExtras *this)
 {
   int pnum,dir;
 
@@ -1062,20 +1069,20 @@ int pugh_SetupBoundingBox(int dim,
       
       for(dir = 0 ; dir < dim; dir++)
       {
-        if (pos[dir] == 0)
+	if (pos[dir] == 0)
 	{
 	  this->lb[pnum][dir] = 0;
 	}
-	else            
+	else		
 	{
 	  this->lb[pnum][dir] = pos[dir]*step[dir] +1 - nghosts[dir];
-	  if(staggertype == PUGH_STAGGER)
-          {
-            this->lb[pnum][dir] --;
-          }
+	  if(stagger == PUGH_STAGGER)
+	  {
+	    this->lb[pnum][dir] --;
+	  }
 	}
 
-        if (pos[dir] == nprocs[dir]-1) 
+	if (pos[dir] == nprocs[dir]-1) 
 	{
 	  this->ub[pnum][dir] = this->nsize[dir]-1;
 	}
@@ -1097,7 +1104,7 @@ int pugh_SetupBoundingBox(int dim,
 
  /*@@
    @routine    pugh_SetupRemoteSizes
-   @date       Mon Nov  8 09:07:27 1999
+   @date       Mon Nov	8 09:07:27 1999
    @author     Tom Goodale
    @desc 
    Determines info about the sizes on each processor.
@@ -1110,13 +1117,13 @@ int pugh_SetupBoundingBox(int dim,
 
 @@*/
 int pugh_SetupRemoteSizes(int dim, 
-                          int periodic,
-                          int staggertype,
-                          int *sh,
-                          int *nghosts,
-                          int total_procs,
-                          int *nprocs,
-                          pGExtras *this)
+			  int periodic,
+			  int stagger,
+			  int *sh,
+			  int *nghosts,
+			  int total_procs,
+			  int *nprocs,
+			  pGExtras *this)
 {
   int pnum;
   int dir;
@@ -1130,9 +1137,301 @@ int pugh_SetupRemoteSizes(int dim,
     for (dir=0;dir<dim;dir++)
     {
       this->rnsize[pnum][dir] = (this->ub[pnum][dir]-this->lb[pnum][dir]+1);
-      this->rnpoints[pnum]  *=  this->rnsize[pnum][dir];
+      this->rnpoints[pnum]  *=	this->rnsize[pnum][dir];
     }
   }
 
   return 0;
 }
+
+ /*@@
+   @routine    SetupPGA
+   @date       Mon Nov  8 16:29:34 1999
+   @author     Tom Goodale
+   @desc 
+   Sets up a new pGA.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+pGA *SetupPGA(void *parent,
+	      pGExtras *extras,
+	      pConnectivity *connectivity,
+	      const char *name,
+	      int id,
+	      int varsize, 
+	      int vtype, 
+	      int stagger) 
+{
+  pGA *this;
+  int i;
+  this = (pGA *)malloc(sizeof(pGA *));
+
+  if(this)
+  {
+    this->extras       = extras;
+    this->connectivity = connectivity;
+    this->parent       = parent;
+    this->varsize      = varsize;
+    this->vtype	       = vtype;
+    this->stagger      = stagger;
+    this->id	       = id;
+
+    this->storage      = PUGH_NOSTORAGE;
+    this->commflag     = PUGH_NOCOMM;
+
+    this->name	      = (char *) malloc((strlen(name)+1)*sizeof(char));
+    this->buffer_sz   = (int *)	 malloc(2*extras->dim*sizeof(int));
+    this->send_buffer = (void **)malloc(2*extras->dim*sizeof(void *));
+    this->recv_buffer = (void **)malloc(2*extras->dim*sizeof(void *));
+    this->docomm      = (int *)	 malloc(2*extras->dim*sizeof(int));
+    this->padddata    = (void *) calloc(1, varsize);
+    this->data	      = this->padddata;
+    
+#ifdef MPI
+    this->sreq = (MPI_Request *)malloc(2*extras->dim*sizeof(MPI_Request));
+    this->rreq = (MPI_Request *)malloc(2*extras->dim*sizeof(MPI_Request));
+#endif
+
+    if(this->name	 &&
+       this->buffer_sz	 &&
+       this->send_buffer &&
+       this->recv_buffer &&
+       this->docomm	 &&
+#ifdef MPI
+       this->ssreq       &&
+       this->rreq        &&
+#endif
+       this->padddata)
+    {
+      strcpy(this->name, name);
+
+      for(i = 0; i < 2*extras->dim; i++)
+      {
+	this->buffer_sz[i]   = 0;
+	this->send_buffer[i] = NULL;
+	this->recv_buffer[i] = NULL;
+	this->docomm[i]	     = 0;
+#ifdef MPI
+        /* Null my send and recieve requests */
+        this->sreq[i] = MPI_REQUEST_NULL;
+        this->rreq[i] = MPI_REQUEST_NULL;
+#endif
+      }
+    }
+    else
+    {
+      free(this->name);
+      free(this->buffer_sz);
+      free(this->send_buffer);
+      free(this->recv_buffer);
+      free(this->docomm);
+      free(this->padddata);
+      
+#ifdef MPI
+      free(this->sreq);
+      free(this->rreq);
+#endif
+      free(this);
+      this = NULL;
+    }
+  }
+
+  return this;
+}
+
+int pugh_EnablePGAStorage(pGA *GA, 
+                          int this_proc,
+                          int zero_memory, 
+                          int padding_active,
+                          int padding_cacheline_bits,
+                          int padding_size,
+                          int padding_address_spacing)
+{
+  int retval;
+  int i;
+  int dir;
+  int sz;
+  int ldir;
+  int special_pad, cache_size;
+  long start;
+
+  retval = 0;
+  if(GA->storage == PUGH_NOSTORAGE)
+  {
+    /* Allocate memory for communication buffers */
+    for(i = 0; i < GA->extras->dim; i++)
+    {
+      dir = i/2;
+
+      sz = GA->extras->nghostzones[dir];      
+      for(ldir = 0; ldir < GA->extras->dim; ldir++)
+      {
+        sz *= GA->extras->lnsize[ldir];
+      }
+      
+      sz /= GA->extras->lnsize[dir];
+
+      if(GA->connectivity->neighbours[this_proc][i] >= 0)
+      {
+        GA->buffer_sz[i]   = sz;
+        GA->send_buffer[i] = malloc(sz * GA->varsize);
+        GA->recv_buffer[i] = malloc(sz * GA->varsize);
+        
+        if(!(GA->send_buffer[i] 
+             && GA->recv_buffer[i]))
+        {
+          retval = -1;
+          break;
+        }
+      }
+      else
+      {
+        GA->buffer_sz[i]   = 0;
+        GA->send_buffer[i] = NULL;
+        GA->recv_buffer[i] = NULL;
+      }
+    }
+
+    if(retval)
+    {
+      for(; i >=0 ; i--)
+      {
+        free(GA->send_buffer[i]);
+        free(GA->recv_buffer[i]);
+        GA->buffer_sz[i] = 0;
+      }
+    }
+
+    /* Now assign memory for the variable itself */
+    if(!retval)
+    {
+      if(GA->padddata)
+      {
+        free(GA->padddata);
+      }
+
+      if(!padding_active)
+      {
+        /* Easy case. */
+        if(zero_memory)
+        {
+          GA->padddata = calloc(GA->extras->npoints, GA->varsize);
+        }
+        else
+        {
+          GA->padddata = malloc(GA->extras->npoints * GA->varsize);
+        }
+
+        GA->data = GA->padddata;
+      }
+      else
+      {
+        if(zero_memory)
+        {
+          GA->padddata = calloc(GA->extras->npoints + padding_size, GA->varsize);
+        }
+        else
+        {
+          GA->padddata = malloc((GA->extras->npoints + padding_size) * GA->varsize);
+        }
+
+           
+        /* Align the actual starting address on a cache line.
+           More specifically align on a unique cache line for
+           each field. 
+        */
+        cache_size = 2 << (padding_cacheline_bits - 1);
+        start = ((long) (GA->padddata) / GA->varsize)%cache_size;
+
+        special_pad = ((GA->id) * padding_address_spacing + cache_size - start)%cache_size;
+
+        GA->data = (char *) GA->padddata + special_pad*GA->varsize;
+
+        if (special_pad > padding_size)
+        {
+          char msg [100];
+          
+          sprintf (msg, "FATAL ERROR: padding size not large enough for cache type "
+                   "specified %d %d %d\n", special_pad, padding_size, cache_size);
+          CCTK_WARN (0, msg);
+          free(msg);
+        }
+
+        retval = 0;
+      }
+    }
+  }
+
+  if (!GA->padddata) 
+  {
+    char *msg = (char *) malloc (80 + strlen (GA->name));
+
+    sprintf (msg, "FATAL ERROR: Cannot allocate data for %s [%d]\n",
+                 GA->name, GA->id);
+    CCTK_WARN (0, msg);
+    free (msg);
+  }
+
+  GA->storage = PUGH_STORAGE;
+
+  return retval;
+
+}
+
+int pugh_DisableGADataStorage(pGA *GA) 
+{
+  int retval;
+  int i;
+
+  if (GA->storage == PUGH_STORAGE) 
+  {
+    if (GA->padddata) 
+    {
+      free(GA->padddata);
+      GA->padddata = NULL;
+      GA->data = NULL;
+    }
+
+    for (i = 0; i < 2*GA->extras->dim; i++)
+    {
+      if (GA->send_buffer[i]) 
+      {
+#ifdef MPI
+        if(GA->sreq[i] != MPI_REQUEST_NULL)
+        {
+          CACTUS_MPI_ERROR(MPI_Request_free(&(GA->sreq[i])));
+        }
+#endif
+        free(GA->send_buffer[i]);
+        GA->send_buffer[i] = NULL;
+      }
+    }
+
+    for (i = 0; i < 2*GA->extras->dim; i++)
+    {
+      if (GA->recv_buffer[i]) 
+      {
+        free(GA->recv_buffer[i]);
+        GA->recv_buffer[i] = NULL;
+      }
+    }
+    
+    GA->padddata = calloc(1, GA->varsize);
+    GA->data = GA->padddata;
+    
+    GA->storage = PUGH_NOSTORAGE;
+    
+    retval = 0;
+  }
+  else
+  {
+    retval = 1;
+  }
+
+  return retval;
+}