Registered reduction routines return an error code, and split reduction.c

into several files


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

1 files changed, 650 insertions, 0 deletions

diff --git a/src/ReductionNorm1.c b/src/ReductionNorm1.c
new file mode 100644
index 0000000..94a9d3f
--- /dev/null
+++ b/src/ReductionNorm1.c
@@ -0,0 +1,650 @@
+ /*@@
+   @file      Reduction_nm1.c
+   @date      Thu Apr  3 11:54:53 1997
+   @author    Thomas Radke, Paul Walker
+   @desc 
+   Various MPI reduction operators. 
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "cctk.h"
+#include "pugh.h"
+
+static char *rcsid = "$Id$";
+
+CCTK_FILEVERSION(CactusPUGH_PUGH_ReductionNorm1_c)
+
+/* some useful macros */
+#ifdef MIN
+#undef MIN
+#endif
+#ifdef ABS
+#undef ABS
+#endif
+#ifdef SQR
+#undef SQR
+#endif
+#define MIN(x, y)   ((x) < (y) ? (x) : (y))
+#define ABS(x)      ((x) < 0   ?-(x) : (x))
+#define SQR(x)      ((x) * (x))
+
+
+/* local function prototypes */
+int PUGH_Norm1_GF  (cGH *GH, int index, int proc, CCTK_REAL *outVal);
+int PUGH_Norm1_Array (cGH *GH, int proc, int nOutVals, int nPoints,
+                      void *inArray, int inType, CCTK_REAL outVals []);
+void copy_real_to_outtype (CCTK_REAL *from, void *to, int to_type, int n);
+int getSizeofElement (int varType);
+
+
+/*@@
+  @routine PUGH_ReduceArrayNorm1
+  @author  Thomas Radke
+  @date    19 Aug 1999
+  @desc
+           Registered PUGH reduction routine for computing the "norm1"
+           of a set of arrays.
+           The arrays are described by their dimensions and variable type.
+           If the number of requested output values equals the number of
+           elements in an input array, global reduction is done element-wise.
+           Otherwise nOutVals must be 1, and global reduction is done on the
+           results of the local reductions.
+           Type casting of the result is provided by specifying the
+           requested output datatype. The intermediate reduction value
+           is always computed as a CCTK_REAL value internally.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nDims
+  @vdesc   number of dimensions of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     dims
+  @vdesc   dimensions of input arrays
+  @vtype   int *
+  @vio     in
+  @endvar 
+  @var     nArrays
+  @vdesc   number of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     inArrays
+  @vdesc   field of input arrays
+  @vtype   void **
+  @vio     in
+  @endvar 
+  @var     inType
+  @vdesc   (common) variable type of input arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nOutVals
+  @vdesc   number of values per output array
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outVals
+  @vdesc   pointer to buffer holding the output values
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     outType
+  @vdesc   (common) variable type of output arrays
+  @vtype   int
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReduceArrayNorm1 (cGH *GH, int proc, int nDims, int dims [],
+                            int nArrays, void *inArrays [], int inType,
+                            int nOutVals, void *outVals, int outType)
+{
+  int retval = 1;
+  int i, elemSize, nPoints;
+  CCTK_REAL *buffer;
+
+
+  nPoints = dims [0];
+  for (i = 1; i < nDims; i++)
+  {
+    nPoints *= dims [i];
+  }
+
+  if (nOutVals != 1) 
+  {
+    if (nOutVals != nPoints) 
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "PUGH_ReduceArrayNorm1: Don't know how to reduce "
+                  "a %d-dimensional array with "
+                  "%d elements to an output array of %d elements",
+                  nDims, nPoints, nOutVals);
+      return -1;
+    }
+    nPoints = 1;
+  }
+
+  if ((elemSize = getSizeofElement (outType)) < 0) 
+  {
+    CCTK_WARN (1, "PUGH_ReduceArrayNorm1: Unsupported return type");
+    return -1;
+  }
+
+  buffer = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+  for (i = 0; i < nArrays; i++) 
+  {
+
+    /* do the reduction on input array [i] */
+    if (PUGH_Norm1_Array (GH, proc, nOutVals, nPoints, inArrays [i], inType,
+                           buffer) == 0) 
+    {
+
+      /* type-cast the result to the requested datatype */
+      if (proc < 0 || proc == CCTK_MyProc (GH))
+      {
+        copy_real_to_outtype (buffer, (char *) outVals + i*nOutVals*elemSize,
+                              outType, nOutVals);
+	retval = 0;
+      }
+    }
+    else
+    {
+      retval = -1;
+    }
+  }
+
+  free (buffer);
+  return retval;
+}
+
+
+
+
+/*@@
+  @routine PUGH_ReduceNorm1
+  @author  Thomas Radke
+  @date    23 Apr 1999
+  @desc
+           PUGH reduction routine to be registered for calculating the norm1.
+           It just goes through the list of variables and calls the appropriate
+           grouptype reduction routine.
+           Global reduction is always done on the results of the local
+           reductions, so nOutVals must be 1.
+  @enddesc
+  @history
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK grid hierarchy
+  @vtype   cGH
+  @vio     in
+  @endvar 
+  @var     proc
+  @vdesc   processor that should receive the result of operation
+           (negative value means all processors receive the result)
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     nOutVars
+  @vdesc   number of requested output elements
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outType
+  @vdesc   type of return value
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     outVals
+  @vdesc   array for the result values to be stored
+  @vtype   void *
+  @vio     in
+  @endvar 
+  @var     numinvars
+  @vdesc   number of variables in the list
+  @vtype   int
+  @vio     in
+  @endvar 
+  @var     varlist
+  @vdesc   list of variables (given as their global indices)
+  @vtype   int *
+  @vio     in
+  @endvar 
+@@*/
+int PUGH_ReduceNorm1 (cGH *GH, int proc,  
+		       int nOutVals,
+                       int outType, 
+		       void *outVals, 
+		       int numinvars, 
+		       int *varlist)
+{
+  int retval=1;
+  int i, elemSize;
+  pGH *pughGH;
+  CCTK_REAL intermediate_result;
+
+  if (nOutVals != 1) 
+  {
+    CCTK_WARN (1, "PUGH_ReduceNorm1: Only one output value per "
+                  "input array allowed");
+    return -2;
+  }
+
+  if ((elemSize = getSizeofElement (outType)) < 0)
+  {
+      CCTK_WARN (1, "PUGH_ReduceNorm1: Unsupported return type");
+      return -3;
+  }
+
+  pughGH = PUGH_pGH (GH);
+
+  for (i = 0; i < numinvars; i++) 
+  {
+    int result;
+
+    result = 1;
+
+    if (CCTK_GroupDimFromVarI(varlist[i])<3)
+    {
+      CCTK_VWarn(1,__LINE__,__FILE__,"PUGH",
+		 "PUGH_ReduceNorm1: Reducing dimension %d not yet supported",
+		 CCTK_GroupDimFromVarI(varlist[i]));
+    }
+    else
+    {
+		 
+      switch (CCTK_GroupTypeFromVarI (varlist [i]))
+      {
+        case GROUP_GF:
+	  result = PUGH_Norm1_GF (GH, varlist [i], proc, &intermediate_result);
+	  break;
+        case GROUP_ARRAY:
+	  CCTK_WARN (1, "PUGH_ReduceNorm1: Reducing arrays not yet supported");
+	  break;
+        case GROUP_SCALAR:
+	  CCTK_WARN (1, "PUGH_ReduceNorm1: Reducing scalars doesn't make sense");
+	  break;
+        default:
+	  CCTK_WARN (1, "PUGH_ReduceNorm1: Invalid group type");
+	  break;
+      }
+    }
+
+    if (result == 0 && (proc < 0 || proc == pughGH->myproc))
+    {
+      copy_real_to_outtype (&intermediate_result, (char *) outVals +
+                            i*nOutVals*elemSize, outType, nOutVals);
+      retval = 0;
+    }
+    else if (result == 0)
+    {
+      retval = -1;
+    }
+  }
+
+  return retval;
+}
+
+
+
+/*********************************************************************/
+/*                        local functions                            */
+/*********************************************************************/
+
+
+/*@@
+   @routine    PUGH_Norm1_Array
+   @date       Aug 19 1999
+   @author     Thomas Radke
+   @desc       Returns the "norm1" of a distributed array with
+               'nPoints' elements. Global reduction is done element-wise
+               (nOutVals == 1) or on the results of the local reductions.
+
+               "norm1" is defined as $\Sigma |a_i| / np$.
+   @enddesc 
+@@*/
+int PUGH_Norm1_Array (cGH *GH, int proc, int nOutVals, int nPoints,
+                      void *inArray, int inType, CCTK_REAL outVals [])
+{
+  int i, j;
+  pGH *pughGH = PUGH_pGH (GH);
+
+
+  for (i = 0; i < nOutVals; i++) 
+  {
+    outVals [i] = 0.0;
+    for (j = 0; j < nPoints; j++) 
+    {
+      if (inType == CCTK_VARIABLE_CHAR)           /* CCTK_CHAR is unsigned */
+      {
+        outVals [i] += ((CCTK_CHAR *) inArray) [i*nPoints + j];
+      }
+      else if (inType == CCTK_VARIABLE_INT)
+      {
+        outVals [i] += ABS (((CCTK_INT *) inArray) [i*nPoints + j]);
+      }
+#ifdef CCTK_INT2
+      else if (inType == CCTK_VARIABLE_INT2)
+      {
+        outVals [i] += ABS (((CCTK_INT2 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_INT4
+      else if (inType == CCTK_VARIABLE_INT4)
+      {
+        outVals [i] += ABS (((CCTK_INT4 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_INT8
+      else if (inType == CCTK_VARIABLE_INT8)
+      {
+        outVals [i] += ABS (((CCTK_INT8 *) inArray) [i*nPoints + j]);
+      }
+#endif
+      else if (inType == CCTK_VARIABLE_REAL)
+      {
+        outVals [i] += ABS (((CCTK_REAL *) inArray) [i*nPoints + j]);
+      }
+#ifdef CCTK_REAL4
+      else if (inType == CCTK_VARIABLE_REAL4)
+      {
+        outVals [i] += ABS (((CCTK_REAL4 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_REAL8
+      else if (inType == CCTK_VARIABLE_REAL8)
+      {
+        outVals [i] += ABS (((CCTK_REAL8 *) inArray) [i*nPoints + j]);
+      }
+#endif
+#ifdef CCTK_REAL16
+      else if (inType == CCTK_VARIABLE_REAL16)
+      {
+        outVals [i] += ABS (((CCTK_REAL16 *) inArray) [i*nPoints + j]);
+      }
+#endif
+      else if (inType == CCTK_VARIABLE_COMPLEX) 
+      {
+        outVals [i] += 
+          sqrt ( SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Re)+
+                 SQR(((CCTK_COMPLEX *) inArray) [i*nPoints + j].Im) );
+      } 
+      else 
+      {
+        CCTK_WARN (1, "PUGH_Norm1_Array: Unknown variable type");
+        return (1);
+      }
+    }
+  }
+
+#ifdef CCTK_MPI
+
+  if (pughGH->nprocs > 1) 
+  {
+    CCTK_REAL *localOutVals;
+
+    localOutVals = (CCTK_REAL *) malloc (nOutVals * sizeof (CCTK_REAL));
+
+    /* outVals[] contains now the local sum */
+    memcpy (localOutVals, outVals, nOutVals * sizeof (CCTK_REAL));
+    if (proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allreduce (localOutVals, outVals, nOutVals,
+                        PUGH_MPI_REAL, MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Reduce (localOutVals, outVals, nOutVals,
+                        PUGH_MPI_REAL, MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    }
+    free (localOutVals);
+
+    nPoints *= pughGH->nprocs;
+  }
+
+#endif
+
+  if (proc < 0 || proc == pughGH->myproc)
+  {
+    for (i = 0; i < nOutVals; i++)
+    {
+      outVals [i] /= nPoints;
+    }
+  }
+
+  return (0);
+}
+
+
+
+/*@@
+   @routine    PUGH_Norm1_GF
+   @date       Thu Apr  3 11:56:05 1997
+   @author     Paul Walker
+   @desc 
+   Returns the average of a distributed grid function.
+   @enddesc 
+@@*/
+int PUGH_Norm1_GF (cGH *GH, int index, int proc, CCTK_REAL *outVal)
+{
+  int timelevel;
+  int ii;
+  CCTK_INT tnp;
+  pGH *pughGH;
+  pGA *GA;
+#ifdef CCTK_MPI
+  int jj, kk;
+  CCTK_INT np;
+  CCTK_REAL local_result;
+#endif
+
+  timelevel = CCTK_NumTimeLevelsFromVarI (index) - 1;
+  if (timelevel > 0)
+  {
+    timelevel--;
+  }
+
+  pughGH = PUGH_pGH (GH);
+  GA = ((pGA ***) pughGH->variables) [index][timelevel];
+
+  *outVal = 0.0;
+
+#ifndef CCTK_MPI
+  for (ii = 0; ii < GA->extras->npoints; ii++)
+  {
+    if (GA->vtype == CCTK_VARIABLE_CHAR)
+    {
+      *outVal += ((CCTK_CHAR *) GA->data) [ii];     /* CCTK_CHAR is unsigned */
+    }
+    else if (GA->vtype == CCTK_VARIABLE_INT)
+    {
+      *outVal += ABS (((CCTK_INT *) GA->data) [ii]);
+    }
+#ifdef CCTK_INT2
+    else if (GA->vtype == CCTK_VARIABLE_INT2)
+    {
+      *outVal += ABS (((CCTK_INT2 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_INT4
+    else if (GA->vtype == CCTK_VARIABLE_INT4)
+    {
+      *outVal += ABS (((CCTK_INT4 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_INT8
+    else if (GA->vtype == CCTK_VARIABLE_INT8)
+    {
+      *outVal += ABS (((CCTK_INT8 *) GA->data) [ii]);
+    }
+#endif
+    else if (GA->vtype == CCTK_VARIABLE_REAL)
+    {
+      *outVal += ABS (((CCTK_REAL *) GA->data) [ii]);
+    }
+#ifdef CCTK_REAL4
+    else if (GA->vtype == CCTK_VARIABLE_REAL4)
+    {
+      *outVal += ABS (((CCTK_REAL4 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_REAL8
+    else if (GA->vtype == CCTK_VARIABLE_REAL8)
+    {
+      *outVal += ABS (((CCTK_REAL8 *) GA->data) [ii]);
+    }
+#endif
+#ifdef CCTK_REAL16
+    else if (GA->vtype == CCTK_VARIABLE_REAL16)
+    {
+      *outVal += ABS (((CCTK_REAL16 *) GA->data) [ii]);
+    }
+#endif
+    else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
+    {
+      *outVal += sqrt( SQR((((CCTK_COMPLEX *) GA->data) [ii]).Re)+
+                       SQR((((CCTK_COMPLEX *) GA->data) [ii]).Im) );
+    } 
+    else 
+    {
+      CCTK_WARN (1, "PUGH_Norm1_GA: Unknown variable type");
+      return (1);
+    }
+  }
+
+  tnp = GA->extras->npoints;
+
+#else
+
+  /* MPI version depends heavily on the ownership array */
+  np  = 0;
+  local_result = 0.0;
+
+  for (kk=GA->extras->ownership[GA->stagger][0][2]; 
+       kk < GA->extras->ownership[GA->stagger][1][2]; 
+       kk ++)
+  {
+    for (jj = GA->extras->ownership[GA->stagger][0][1]; 
+         jj < GA->extras->ownership[GA->stagger][1][1]; 
+         jj ++)
+    {
+      for (ii = GA->extras->ownership[GA->stagger][0][0]; 
+           ii < GA->extras->ownership[GA->stagger][1][0]; 
+           ii++, np++)
+      {
+        if (GA->vtype == CCTK_VARIABLE_CHAR)
+        {
+          local_result +=              /* CCTK_CHAR is unsigned */
+	    ((CCTK_CHAR *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)];
+        }
+        else if (GA->vtype == CCTK_VARIABLE_INT)
+        {
+          local_result +=
+	    ABS (((CCTK_INT *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#ifdef CCTK_INT2
+        else if (GA->vtype == CCTK_VARIABLE_INT2)
+        {
+          local_result +=
+	    ABS (((CCTK_INT2 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_INT4
+        else if (GA->vtype == CCTK_VARIABLE_INT4)
+        {
+          local_result +=
+	    ABS (((CCTK_INT4 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_INT8
+        else if (GA->vtype == CCTK_VARIABLE_INT8)
+        {
+          local_result +=
+	    ABS (((CCTK_INT8 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+        else if (GA->vtype == CCTK_VARIABLE_REAL)
+        {
+          local_result +=
+	    ABS (((CCTK_REAL *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#ifdef CCTK_REAL4
+        else if (GA->vtype == CCTK_VARIABLE_REAL4)
+        {
+          local_result +=
+	    ABS (((CCTK_REAL4 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_REAL8
+        else if (GA->vtype == CCTK_VARIABLE_REAL8)
+        {
+          local_result +=
+	    ABS (((CCTK_REAL8 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+#ifdef CCTK_REAL16
+        else if (GA->vtype == CCTK_VARIABLE_REAL16)
+        {
+          local_result +=
+	    ABS (((CCTK_REAL16 *) GA->data) [DATINDEX (GA->extras, ii, jj, kk)]);
+        }
+#endif
+        else if (GA->vtype == CCTK_VARIABLE_COMPLEX) 
+        {
+          local_result +=
+            sqrt( SQR((((CCTK_COMPLEX *)GA->data)
+		       [DATINDEX (GA->extras,ii,jj,kk)]).Re)+
+		  SQR((((CCTK_COMPLEX *)GA->data)
+		       [DATINDEX (GA->extras,ii,jj,kk)]).Im) );  
+        } 
+        else 
+        {
+          CCTK_WARN (1, "PUGH_Norm1_GF: Unknown variable type");
+          return (1);
+        }
+      }
+    }
+  }
+
+  if (proc < 0) 
+  {
+    CACTUS_MPI_ERROR (MPI_Allreduce (&local_result, outVal, 1, PUGH_MPI_REAL,
+                      MPI_SUM, pughGH->PUGH_COMM_WORLD));
+    CACTUS_MPI_ERROR (MPI_Allreduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
+                      pughGH->PUGH_COMM_WORLD));
+  } 
+  else 
+  {
+    CACTUS_MPI_ERROR (MPI_Reduce (&local_result, outVal, 1, PUGH_MPI_REAL,
+                      MPI_SUM, proc, pughGH->PUGH_COMM_WORLD));
+    CACTUS_MPI_ERROR (MPI_Reduce (&np, &tnp, 1, PUGH_MPI_INT, MPI_SUM,
+                      proc, pughGH->PUGH_COMM_WORLD));
+  }
+#endif
+ 
+  if (proc < 0 || proc == pughGH->myproc)
+  {
+    *outVal /= tnp;
+  }
+
+  return (0);
+}
+
+