This commit was generated by cvs2svn to compensate for changes in r2, which

included commits to RCS files with non-trunk default branches.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHReduce/trunk@3 d60812e6-3970-4df4-986e-c251b06effeb

1 files changed, 461 insertions, 0 deletions

diff --git a/src/Reduction.c b/src/Reduction.c
new file mode 100644
index 0000000..80266d8
--- /dev/null
+++ b/src/Reduction.c
@@ -0,0 +1,461 @@
+ /*@@
+   @file      Reduction.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc 
+   Various MPI reduction operators. 
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include <stdlib.h>
+
+#include "pugh_reductions.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGHReduce_Reduction_c)
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+static int PUGH_ReductionGA (cGH *GH, 
+			     int index, 
+			     int proc, 
+			     CCTK_REAL *outval,
+                             reduction_fn_t reduction_fn);
+
+static int copy_real_to_outtype (int num_elems,
+                                 CCTK_REAL inarray [/* num_elems */],
+                                 int outtype,
+                                 void *outarray     /* [num_elems] */);
+
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    PUGH_ReductionArrays
+   @date       
+   @author     
+   @desc 
+
+   @enddesc 
+   @calls      
+   @calledby   
+   @history 
+ 
+   @endhistory 
+   @var     
+   @vdesc   
+   @vtype   
+   @vio     
+   @vcomment 
+ 
+   @endvar 
+
+@@*/
+
+int PUGH_ReductionArrays (cGH *GH,
+                          int proc,
+                          int num_dims,
+                          int dims[/* num_dims */],
+                          int intype,
+                          int num_inarrays,
+                          void *inarrays[/* num_inarrays */],
+                          int outtype,
+                          int num_outvals,
+                          void *outvals /* [num_outvals] */,
+                          reduction_fn_t reduction_fn)
+{
+  int retval;
+  int i;
+  int num_points;
+  int from[1];
+  int to[1];
+  int iterator[1];
+  int points_per_dim[1];
+  int *intypes;
+  CCTK_REAL *buffer;
+
+  points_per_dim[0] = 1;
+  from[0] = 0;
+  to[0] = dims[0];
+
+  for (i = 1; i < num_dims; i++)
+  {
+    to[0] *= dims[i];
+  }
+
+  if (num_outvals != 1)
+  {
+    if (num_outvals != to[0])
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "PUGH_ReductionArrays: Don't know how to reduce "
+                  "a %d-dimensional array with %d elements "
+                  "to an output array of %d elements",
+                  num_dims, to[0], num_outvals);
+      return (-1);
+    }
+    to[0] = 1;
+  }
+  num_points = to[0] * CCTK_nProcs (GH);
+
+  /* set the array types to intype */
+  /* FIXME: could allow to pass in arrays of different types now !!! */
+  intypes = (int *) malloc (num_inarrays * sizeof (int));
+  for (i = 0; i < num_inarrays; i++)
+  {
+    intypes[i] = intype;
+  }
+
+  buffer = (CCTK_REAL *) malloc (num_outvals * sizeof (CCTK_REAL));
+
+  /* do the reduction on the input arrays */
+  retval = reduction_fn (GH, 
+			 proc, 
+			 num_dims, 
+			 from, 
+			 to, 
+			 iterator, 
+			 points_per_dim,
+                         num_points, 
+			 num_inarrays, 
+			 intypes, 
+			 inarrays,
+                         num_outvals, 
+			 buffer);
+
+  if (retval == 0 && (proc < 0 || proc == CCTK_MyProc (GH)))
+  {
+    /* type-cast the result to the requested datatype */
+    retval = copy_real_to_outtype (num_inarrays * num_outvals, 
+				   buffer,
+                                   outtype, 
+				   outvals);
+  }
+
+  free (intypes);
+  free (buffer);
+
+  return (retval);
+}
+
+
+ /*@@
+   @routine    PUGH_ReductionArrays
+   @date       
+   @author     
+   @desc 
+
+   @enddesc 
+   @calls      
+   @calledby   
+   @history 
+ 
+   @endhistory 
+   @var     
+   @vdesc   
+   @vtype   
+   @vio     
+   @vcomment 
+ 
+   @endvar 
+
+@@*/
+
+int PUGH_ReductionGVs (cGH *GH,
+                       int proc,
+                       int num_invars,
+                       int invars[/* num_invars */],
+                       int outtype,
+                       int num_outvals,
+                       void *outvals /* [num_outvals] */,
+                       reduction_fn_t reduction_fn)
+{
+  int i, myproc, outtypesize, this_retval, retval;
+  CCTK_REAL result;
+
+
+  if (num_outvals != 1)
+  {
+    CCTK_WARN (1, "PUGH_ReductionGVs: Only one output value per "
+                  "input array allowed");
+    return (-1);
+  }
+
+  outtypesize = CCTK_VarTypeSize (outtype);
+  if (outtypesize <= 0)
+  {
+    CCTK_WARN (1, "PUGH_ReductionGVs: Invalid output variable type");
+    return (-1);
+  }
+
+  retval = 0;
+  myproc = CCTK_MyProc (GH);
+
+  for (i = 0; i < num_invars; i++)
+  {
+    switch (CCTK_GroupTypeFromVarI (invars[i]))
+    {
+      case CCTK_GF:
+      case CCTK_ARRAY:
+        this_retval = PUGH_ReductionGA (GH, invars[i], proc, &result,
+                                        reduction_fn);
+        break;
+
+      case CCTK_SCALAR:
+        CCTK_WARN (1, "PUGH_ReductionGVs: Reducing scalars doesn't make sense");
+        this_retval = -1;
+        break;
+
+      default:
+        CCTK_WARN (1, "PUGH_ReductionGVs: Invalid variable index");
+        this_retval = -1;
+        break;
+    }
+
+    if (this_retval == 0 && (proc < 0 || proc == myproc))
+    {
+      /* type-cast the result to the requested datatype */
+      this_retval = copy_real_to_outtype (1, &result, outtype,
+                                          (char *) outvals + i*outtypesize);
+    }
+
+    retval |= this_retval;
+  }
+
+  return (retval);
+}
+
+/********************************************************************
+ *********************     Local Routines   *************************
+ ********************************************************************/
+
+ /*@@
+   @routine    PUGH_ReductionGA
+   @date       
+   @author     
+   @desc 
+
+   @enddesc 
+   @calls      
+   @calledby   
+   @history 
+ 
+   @endhistory 
+   @var     
+   @vdesc   
+   @vtype   
+   @vio     
+   @vcomment 
+ 
+   @endvar 
+
+@@*/
+
+static int PUGH_ReductionGA (cGH *GH, int index, int proc, CCTK_REAL *outval,
+                             reduction_fn_t reduction_fn)
+{
+  int i, retval, stagger_index, num_points;
+  pGA *GA;
+  int *from, *to, *iterator, *points_per_dim;
+
+
+  /* get the GA structure for the current timelevel */
+  i = CCTK_NumTimeLevelsFromVarI (index) - 1;
+  if (i > 0)
+  {
+    i--;
+  }
+  GA = ((pGA ***) PUGH_pGH (GH)->variables)[index][i];
+
+  /* allocate temporary vectors */
+  from           = (int *) malloc (4 * GA->connectivity->dim * sizeof (int));
+  to             = from + 1*GA->connectivity->dim;
+  iterator       = from + 2*GA->connectivity->dim;
+  points_per_dim = from + 3*GA->connectivity->dim;
+
+  /* get the start and endpoint to iterate over, the total number of points,
+     and the points_per_dim[] vector */
+  num_points = 1;
+  points_per_dim[0] = 1;
+  for (i = 0; i < GA->connectivity->dim; i++)
+  {
+    stagger_index = CCTK_StaggerDirIndex (i, GA->stagger);
+
+    from[i] = GA->extras->ownership[stagger_index][0][i];
+    to[i]   = GA->extras->ownership[stagger_index][1][i];
+
+    if (stagger_index == 0)
+    {
+      num_points *= GA->extras->nsize[i];
+    }
+    else
+    {
+      num_points *= GA->extras->nsize[i] - 1;
+    }
+
+    if (i > 0)
+    {
+      points_per_dim[i] = points_per_dim[i-1] * GA->extras->lnsize[i-1];
+    }
+  }
+
+  /* now do the reduction */
+  retval = reduction_fn (GH, proc, GA->connectivity->dim, from, to, iterator,
+                         points_per_dim, num_points, 1, &GA->vtype, &GA->data,
+                         1, outval);
+
+  /* free temporary vectors */
+  free (from);
+
+  return (retval);
+}
+
+
+ /*@@
+   @routine    copy_real_to_outtype
+   @date       
+   @author     
+   @desc 
+
+   @enddesc 
+   @calls      
+   @calledby   
+   @history 
+ 
+   @endhistory 
+   @var     
+   @vdesc   
+   @vtype   
+   @vio     
+   @vcomment 
+ 
+   @endvar 
+
+@@*/
+
+static int copy_real_to_outtype (int num_elems,
+                                 CCTK_REAL inarray[/* num_elems */],
+                                 int outtype,
+                                 void *outarray /* [num_elems] */)
+{
+  int i, retval;
+
+
+  retval = 0;
+
+  if (outtype == CCTK_VARIABLE_CHAR)
+  {
+    CCTK_CHAR *_outarray = (CCTK_CHAR *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_CHAR) inarray[i];
+    }
+  }
+  else if (outtype == CCTK_VARIABLE_INT)
+  {
+    CCTK_INT *_outarray = (CCTK_INT *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_INT) inarray[i];
+    }
+  }
+#ifdef CCTK_INT2
+  else if (outtype == CCTK_VARIABLE_INT2)
+  {
+    CCTK_INT2 *_outarray = (CCTK_INT2 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_INT2) inarray[i];
+    }
+  }
+#endif
+#ifdef CCTK_INT4
+  else if (outtype == CCTK_VARIABLE_INT4)
+  {
+    CCTK_INT4 *_outarray = (CCTK_INT4 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_INT4) inarray[i];
+    }
+  }
+#endif
+#ifdef CCTK_INT8
+  else if (outtype == CCTK_VARIABLE_INT8)
+  {
+    CCTK_INT8 *_outarray = (CCTK_INT8 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_INT8) inarray[i];
+    }
+  }
+#endif
+  else if (outtype == CCTK_VARIABLE_REAL)
+  {
+    CCTK_REAL *_outarray = (CCTK_REAL *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_REAL) inarray[i];
+    }
+  }
+#ifdef CCTK_REAL4
+  else if (outtype == CCTK_VARIABLE_REAL4)
+  {
+    CCTK_REAL4 *_outarray = (CCTK_REAL4 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_REAL4) inarray[i];
+    }
+  }
+#endif
+#ifdef CCTK_REAL8
+  else if (outtype == CCTK_VARIABLE_REAL8)
+  {
+	CCTK_REAL8 *_outarray = (CCTK_REAL8 *) outarray;
+
+
+	for (i = 0; i < num_elems; i++)
+	{
+	  _outarray[i] = (CCTK_REAL8) inarray[i];
+	}
+  }
+#endif
+#ifdef CCTK_REAL16
+  else if (outtype == CCTK_VARIABLE_REAL16)
+  {
+    CCTK_REAL16 *_outarray = (CCTK_REAL16 *) outarray;
+
+
+    for (i = 0; i < num_elems; i++)
+    {
+      _outarray[i] = (CCTK_REAL16) inarray[i];
+    }
+  }
+#endif
+  else
+  {
+    CCTK_WARN (1, "copy_real_to_outtype: Unsupported output type");
+    retval = -1;
+  }
+
+  return (retval);
+}