Added initial stab at reduction operators for Carpet.

darcs-hash:20011211120856-07bb3-0f661df0de1c74d8064bf095db9f5abd93efe9a4.gz

8 files changed, 259 insertions, 850 deletions

diff --git a/Carpet/CarpetReduce/README b/Carpet/CarpetReduce/README
index 8692a2930..e5c64331a 100644
--- a/Carpet/CarpetReduce/README
+++ b/Carpet/CarpetReduce/README
@@ -1,29 +1,7 @@
 Cactus Code Thorn CarpetReduce
 Authors    : Erik Schnetter <schnetter@uni-tuebingen.de>
-CVS info   : $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/README,v 1.3 2004/06/14 07:01:21 schnetter Exp $
+CVS info   : $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/README,v 1.1 2001/12/11 13:08:56 schnetter Exp $
 --------------------------------------------------------------------------
 
 Purpose of the thorn:
 
-This thorn provides parallel reduction operators for Carpet.
-
-
-
-This thorn now uses a weight function.  This makes it possible to
-perform physically meaningful spatial reduction operations.  The
-weight is 1 for all "normal" grid points.
-
-The weight is set to 0 on symmetry and possible the outer boundary,
-and it might be set to 1/2 on the edge of the boundary.  Setting this
-depends on the coordinate thorn, and currently works only when the
-coordinates are defined via CoordBase.
-
-The weight is also reduced or set to 0 on coarser grids that are
-overlaid by finer grid.
-
-The weight should also be reduced or set to 0 near and in excised
-regions.  This should happen in conjunction with an excision boundary
-thorn.
-
-This weigth function should probably be extracted into its own thorn
-MaskBase, so that many thorns can easily operate on it.
diff --git a/Carpet/CarpetReduce/interface.ccl b/Carpet/CarpetReduce/interface.ccl
index 6c14b6bae..f479b5272 100644
--- a/Carpet/CarpetReduce/interface.ccl
+++ b/Carpet/CarpetReduce/interface.ccl
@@ -1,27 +1,6 @@
 # Interface definition for thorn CarpetReduce
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/interface.ccl,v 1.9 2004/06/21 12:37:07 hawke Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/interface.ccl,v 1.1 2001/12/11 13:08:56 schnetter Exp $
 
-IMPLEMENTS: reduce
+implements: reduce
 
-uses include header: dist.hh
-uses include header: vect.hh
-
-uses include header: carpet.hh
-
-
-
-CCTK_INT FUNCTION \
-    SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST IN cctkGH)
-REQUIRES FUNCTION SymmetryTableHandleForGrid
-
-CCTK_INT FUNCTION GetBoundarySpecification \
-  (CCTK_INT IN size, \
-   CCTK_INT OUT ARRAY nboundaryzones, \
-   CCTK_INT OUT ARRAY is_internal, \
-   CCTK_INT OUT ARRAY is_staggered, \
-   CCTK_INT OUT ARRAY shiftout)
-REQUIRES FUNCTION GetBoundarySpecification
-
-
-
-REAL weight TYPE=gf TAGS='prolongation="none"' "Weight function"
+inherits: CarpetLib driver
diff --git a/Carpet/CarpetReduce/param.ccl b/Carpet/CarpetReduce/param.ccl
index 0842ac0db..70fbf03e5 100644
--- a/Carpet/CarpetReduce/param.ccl
+++ b/Carpet/CarpetReduce/param.ccl
@@ -1,6 +1,2 @@
 # Parameter definitions for thorn CarpetReduce
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/param.ccl,v 1.2 2004/06/14 07:01:21 schnetter Exp $
-
-BOOLEAN verbose "Produce screen output while running"
-{
-} "no"
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/param.ccl,v 1.1 2001/12/11 13:08:56 schnetter Exp $
diff --git a/Carpet/CarpetReduce/schedule.ccl b/Carpet/CarpetReduce/schedule.ccl
index db1530ae7..82626dee2 100644
--- a/Carpet/CarpetReduce/schedule.ccl
+++ b/Carpet/CarpetReduce/schedule.ccl
@@ -1,44 +1,7 @@
 # Schedule definitions for thorn CarpetReduce
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/schedule.ccl,v 1.4 2004/08/02 11:43:35 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/schedule.ccl,v 1.1 2001/12/11 13:08:56 schnetter Exp $
 
 schedule CarpetReduceStartup at STARTUP
 {
   LANG: C
 } "Startup routine"
-
-
-
-# This might move to MaskBase
-STORAGE: weight
-
-SCHEDULE GROUP MaskBase_SetupMask AT basegrid
-{
-} "Set up the weight function"
-
-SCHEDULE GROUP MaskBase_SetupMask AT postregrid
-{
-} "Set up the weight function"
-
-SCHEDULE MaskBase_InitMask IN MaskBase_SetupMask
-{
-  LANG: C
-  OPTIONS: global loop-local
-} "Initialise the weight function"
-
-SCHEDULE GROUP SetupMask IN MaskBase_SetupMask AFTER MaskBase_InitMask
-{
-} "Set up the weight function (schedule other routines in here)"
-
-# This might move to CoordBase
-SCHEDULE CoordBase_SetupMask IN SetupMask
-{
-  LANG: C
-  OPTIONS: global loop-local
-} "Set up the outer boundaries of the weight function"
-
-# This might move to CarpetMask
-SCHEDULE CarpetMaskSetup IN SetupMask
-{
-  LANG: C
-  OPTIONS: global loop-singlemap
-} "Set up the weight function for the restriction regions"
diff --git a/Carpet/CarpetReduce/src/make.code.defn b/Carpet/CarpetReduce/src/make.code.defn
index afc3646db..3738a14d4 100644
--- a/Carpet/CarpetReduce/src/make.code.defn
+++ b/Carpet/CarpetReduce/src/make.code.defn
@@ -1,8 +1,8 @@
 # Main make.code.defn file for thorn CarpetReduce
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/make.code.defn,v 1.2 2004/06/14 07:01:21 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/make.code.defn,v 1.1 2001/12/11 13:08:57 schnetter Exp $
 
 # Source files in this directory
-SRCS = mask_carpet.cc mask_coords.c mask_init.c reduce.cc
+SRCS = reduce.cc
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/Carpet/CarpetReduce/src/reduce.cc b/Carpet/CarpetReduce/src/reduce.cc
index 622d63291..8e1c7821f 100644
--- a/Carpet/CarpetReduce/src/reduce.cc
+++ b/Carpet/CarpetReduce/src/reduce.cc
@@ -1,254 +1,33 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.43 2004/08/02 11:43:35 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.1 2001/12/11 13:08:58 schnetter Exp $
 
 #include <assert.h>
-#include <float.h>
 #include <limits.h>
 #include <math.h>
-#include <stdio.h>
 #include <stdlib.h>
 
-#include <complex>
-#include <limits>
-#include <vector>
-
 #include <mpi.h>
 
 #include "cctk.h"
 
-#include "dist.hh"
-#include "vect.hh"
+#include "Carpet/CarpetLib/src/dist.hh"
 
-#include "carpet.hh"
+#include "Carpet/Carpet/src/carpet.hh"
 
 #include "reduce.hh"
 
-extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.43 2004/08/02 11:43:35 schnetter Exp $";
-  CCTK_FILEVERSION(Carpet_CarpetReduce_reduce_cc);
-}
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.1 2001/12/11 13:08:58 schnetter Exp $";
 
 
 
 namespace CarpetReduce {
   
-  using namespace std;
   using namespace Carpet;
   
   
   
-  // Helper functions and types
-  
-  // The minimum of two values
-  template<typename T> inline T
-  mymin (const T x, const T y)
-  {
-    return min(x, y);
-  }
-  
-  // The maximum of two values
-  template<typename T> inline T
-  mymax (const T x, const T y)
-  {
-    return max(x, y);
-  }
-  
-  // Square root
-  template<typename T> inline T
-  mysqrt (const T x)
-  {
-    return sqrt(x);
-  }
-  
-  
-  
-  // Properties of numeric types
-  template<typename T>
-  struct my_numeric_limits {
-    
-    // The smallest possible value
-    static T min ()
-    {
-      return mymin (numeric_limits<T>::min(), -numeric_limits<T>::max());
-    }
-    
-    // The largest possible value
-    static T max ()
-    {
-      return mymax (numeric_limits<T>::max(), -numeric_limits<T>::min());
-    }
-    
-  };
-  
-  
-  
-  // Provide for each Cactus type a "good" type, translating between
-  // CCTK_COMPLEX* and complex<CCTK_REAL*>
-  template<typename T>
-  struct typeconv {
-    typedef T goodtype;
-    typedef T badtype;
-  };
-  
-  
-  
-  // Overload the above helper functions and types for complex values
-  
-#ifdef CCTK_REAL4
-  
-  template<> inline complex<CCTK_REAL4>
-  mymin (const complex<CCTK_REAL4> x, const complex<CCTK_REAL4> y)
-  {
-    return complex<CCTK_REAL4> (mymin(x.real(), y.real()),
-                                mymin(x.imag(), y.imag()));
-  }
-  
-  template<> inline complex<CCTK_REAL4>
-  mymax (const complex<CCTK_REAL4> x, const complex<CCTK_REAL4> y)
-  {
-    return complex<CCTK_REAL4> (mymax(x.real(), y.real()),
-                                mymax(x.imag(), y.imag()));
-  }
-  
-  template<>
-  struct my_numeric_limits<complex<CCTK_REAL4> > {
-    static complex<CCTK_REAL4> min ()
-    {
-      return complex<CCTK_REAL4> (my_numeric_limits<CCTK_REAL4>::min(),
-                                  my_numeric_limits<CCTK_REAL4>::min());
-    }
-    static complex<CCTK_REAL4> max ()
-    {
-      return complex<CCTK_REAL4> (my_numeric_limits<CCTK_REAL4>::max(),
-                                  my_numeric_limits<CCTK_REAL4>::max());
-    }
-  };
-
-  template<>
-  struct typeconv<CCTK_COMPLEX8> {
-    typedef complex<CCTK_REAL4> goodtype;
-    typedef CCTK_COMPLEX8 badtype;
-  };
-  
-#endif
-  
-#ifdef CCTK_REAL8
-  
-  template<> inline complex<CCTK_REAL8>
-  mymin (const complex<CCTK_REAL8> x, const complex<CCTK_REAL8> y)
-  {
-    return complex<CCTK_REAL8> (mymin(x.real(), y.real()),
-                                mymin(x.imag(), y.imag()));
-  }
-  
-  template<> inline complex<CCTK_REAL8>
-  mymax (const complex<CCTK_REAL8> x, const complex<CCTK_REAL8> y)
-  {
-    return complex<CCTK_REAL8> (mymax(x.real(), y.real()),
-                                mymax(x.imag(), y.imag()));
-  }
-  
-  template<>
-  struct my_numeric_limits<complex<CCTK_REAL8> > {
-    static complex<CCTK_REAL8> min ()
-    {
-      return complex<CCTK_REAL8> (my_numeric_limits<CCTK_REAL8>::min(),
-                                  my_numeric_limits<CCTK_REAL8>::min());
-    }
-    static complex<CCTK_REAL8> max ()
-    {
-      return complex<CCTK_REAL8> (my_numeric_limits<CCTK_REAL8>::max(),
-                                  my_numeric_limits<CCTK_REAL8>::max());
-    }
-  };
-
-  template<>
-  struct typeconv<CCTK_COMPLEX16> {
-    typedef complex<CCTK_REAL8> goodtype;
-    typedef CCTK_COMPLEX16 badtype;
-  };
-  
-#endif
-  
-#ifdef CCTK_REAL16
-  
-  template<> inline complex<CCTK_REAL16>
-  mymin (const complex<CCTK_REAL16> x, const complex<CCTK_REAL16> y)
-  {
-    return complex<CCTK_REAL16> (mymin(x.real(), y.real()),
-                                 mymin(x.imag(), y.imag()));
-  }
-  
-  template<> inline complex<CCTK_REAL16>
-  mymax (const complex<CCTK_REAL16> x, const complex<CCTK_REAL16> y)
-  {
-    return complex<CCTK_REAL16> (mymax(x.real(), y.real()),
-                                 mymax(x.imag(), y.imag()));
-  }
-  
-  template<>
-  struct my_numeric_limits<complex<CCTK_REAL16> > {
-    static complex<CCTK_REAL16> min ()
-    {
-      return complex<CCTK_REAL16> (my_numeric_limits<CCTK_REAL16>::min(),
-                                   my_numeric_limits<CCTK_REAL16>::min());
-    }
-    static complex<CCTK_REAL16> max ()
-    {
-      return complex<CCTK_REAL16> (my_numeric_limits<CCTK_REAL16>::max(),
-                                   my_numeric_limits<CCTK_REAL16>::max());
-    }
-  };
-
-  template<>
-  struct typeconv<CCTK_COMPLEX32> {
-    typedef complex<CCTK_REAL16> goodtype;
-    typedef CCTK_COMPLEX32 badtype;
-  };
-  
-#endif
-  
-  
-  
-  // Provide a square root function for integer values
-  
-#ifdef CCTK_INT1
-  template<> inline CCTK_INT1
-  mysqrt (const CCTK_INT1 x)
-  {
-    return sqrt((CCTK_REAL)x);
-  }
-#endif
-  
-#ifdef CCTK_INT2
-  template<> inline CCTK_INT2
-  mysqrt (const CCTK_INT2 x)
-  {
-    return sqrt((CCTK_REAL)x);
-  }
-#endif
-  
-#ifdef CCTK_INT4
-  template<> inline CCTK_INT4
-  mysqrt (const CCTK_INT4 x)
-  {
-    return sqrt((CCTK_REAL)x);
-  }
-#endif
-  
-#ifdef CCTK_INT8
-  template<> inline CCTK_INT8
-  mysqrt (const CCTK_INT8 x)
-  {
-    return sqrt((CCTK_REAL)x);
-  }
-#endif
-  
-  
-  
   // Poor man's RTTI
-  enum ared { do_count, do_minimum, do_maximum, do_product, do_sum,
-              do_sum_abs, do_sum_squared, do_average, do_norm1, do_norm2,
-              do_norm_inf };
+  enum ared { do_count, do_minimum, do_maximum, do_sum, do_sum_abs,
+	      do_mean, do_norm1, do_norm2, do_norm_inf };
   
   
   
@@ -262,146 +41,150 @@ namespace CarpetReduce {
   
   // count: count the number of grid points
   struct count : reduction {
-    count () { }
     ared thered () const { return do_count; }
     bool uses_cnt () const { return false; }
     template<class T>
     struct op {
-      static inline void initialise (T& accum) { accum = T(0); }
-      static inline void reduce (T& accum, const T& val, const CCTK_REAL weight) { accum += T(weight); }
+      static inline void initialise (T& accum) { accum = 0; }
+      static inline void reduce (T& accum, const T& val) { ++accum; }
       static inline void finalise (T& accum, const T& cnt) { }
     };
     MPI_Op mpi_op () const { return MPI_SUM; }
   };
   
   struct minimum : reduction {
-    minimum () { }
     ared thered () const { return do_minimum; }
     bool uses_cnt () const { return false; }
     template<class T>
     struct op {
-      static inline void initialise (T& accum) { accum = my_numeric_limits<T>::max(); }
-      static inline void reduce (T& accum, const T& val, const CCTK_REAL weight) { if (weight>0) accum = mymin(accum, val); }
+      static inline void initialise (T& accum) { abort(); }
+      static inline void reduce (T& accum, const T& val) { accum = min(accum, val); }
       static inline void finalise (T& accum, const T& cnt) { }
     };
     MPI_Op mpi_op () const { return MPI_MIN; }
   };
+  template<> inline void minimum::op<char>         ::initialise (char         & accum) { accum = CHAR_MAX; }
+  template<> inline void minimum::op<signed char>  ::initialise (signed char  & accum) { accum = SCHAR_MAX; }
+  template<> inline void minimum::op<unsigned char>::initialise (unsigned char& accum) { accum = UCHAR_MAX; }
+  template<> inline void minimum::op<short>        ::initialise (short        & accum) { accum = SHRT_MAX; }
+  template<> inline void minimum::op<int>          ::initialise (int          & accum) { accum = INT_MAX; }
+  template<> inline void minimum::op<long>         ::initialise (long         & accum) { accum = LONG_MAX; }
+#ifdef LLONG_MAX
+  template<> inline void minimum::op<long long>    ::initialise (long long    & accum) { accum = LLONG_MAX; }
+#endif
+#ifdef HUGE_VALF
+  template<> inline void minimum::op<float>        ::initialise (float        & accum) { accum = HUGE_VALF; }
+#endif
+  template<> inline void minimum::op<double>       ::initialise (double       & accum) { accum = HUGE_VAL; }
+#ifdef HUGE_VALL
+  template<> inline void minimum::op<long double>  ::initialise (long double  & accum) { accum = HUGE_VALL; }
+#endif
   
   struct maximum : reduction {
-    maximum () { }
     ared thered () const { return do_maximum; }
     bool uses_cnt () const { return false; }
     template<class T>
     struct op {
-      static inline void initialise (T& accum) { accum = my_numeric_limits<T>::min(); }
-      static inline void reduce (T& accum, const T& val, const CCTK_REAL weight) { if (weight>0) accum = mymax(accum, val); }
+      static inline void initialise (T& accum) { abort(); }
+      static inline void reduce (T& accum, const T& val) { accum = max(accum, val); }
       static inline void finalise (T& accum, const T& cnt) { }
     };
     MPI_Op mpi_op () const { return MPI_MAX; }
   };
-  
-  struct product : reduction {
-    product () { }
-    ared thered () const { return do_product; }
-    bool uses_cnt () const { return false; }
-    template<class T>
-    struct op {
-      static inline void initialise (T& accum) { accum = T(1); }
-      static inline void reduce (T& accum, const T& val, const CCTK_REAL weight) { if (weight>0) accum *= weight==1 ? val : pow(val,weight); }
-      static inline void finalise (T& accum, const T& cnt) { }
-    };
-    MPI_Op mpi_op () const { return MPI_PROD; }
-  };
+  template<> inline void maximum::op<char>         ::initialise (char         & accum) { accum = CHAR_MIN; }
+  template<> inline void maximum::op<signed char>  ::initialise (signed char  & accum) { accum = SCHAR_MIN; }
+  template<> inline void maximum::op<unsigned char>::initialise (unsigned char& accum) { accum = 0; }
+  template<> inline void maximum::op<short>        ::initialise (short        & accum) { accum = SHRT_MIN; }
+  template<> inline void maximum::op<int>          ::initialise (int          & accum) { accum = INT_MIN; }
+  template<> inline void maximum::op<long>         ::initialise (long         & accum) { accum = LONG_MIN; }
+#ifdef LLONG_MIN
+  template<> inline void maximum::op<long long>    ::initialise (long long    & accum) { accum = LLONG_MIN; }
+#endif
+#ifdef HUGE_VALF
+  template<> inline void maximum::op<float>        ::initialise (float        & accum) { accum = -HUGE_VALF; }
+#endif
+  template<> inline void maximum::op<double>       ::initialise (double       & accum) { accum = -HUGE_VAL; }
+#ifdef HUGE_VALL
+  template<> inline void maximum::op<long double>  ::initialise (long double  & accum) { accum = -HUGE_VALL; }
+#endif
   
   struct sum : reduction {
-    sum () { }
     ared thered () const { return do_sum; }
     bool uses_cnt () const { return false; }
     template<class T>
     struct op {
-      static inline void initialise (T& accum) { accum = T(0); }
-      static inline void reduce (T& accum, const T& val, const CCTK_REAL weight) { if (weight>0) accum += weight*val; }
+      static inline void initialise (T& accum) { accum = 0; }
+      static inline void reduce (T& accum, const T& val) { accum += val; }
       static inline void finalise (T& accum, const T& cnt) { }
     };
     MPI_Op mpi_op () const { return MPI_SUM; }
   };
   
   struct sum_abs : reduction {
-    sum_abs () { }
     ared thered () const { return do_sum_abs; }
     bool uses_cnt () const { return false; }
     template<class T>
     struct op {
-      static inline void initialise (T& accum) { accum = T(0); }
-      static inline void reduce (T& accum, const T& val, const CCTK_REAL weight) { if (weight>0) accum += weight*abs(val); }
+      static inline void initialise (T& accum) { accum = 0; }
+      static inline void reduce (T& accum, const T& val) { accum += abs(val); }
       static inline void finalise (T& accum, const T& cnt) { }
     };
     MPI_Op mpi_op () const { return MPI_SUM; }
   };
   
-  struct sum_squared : reduction {
-    sum_squared () { }
-    ared thered () const { return do_sum_squared; }
-    bool uses_cnt () const { return false; }
-    template<class T>
-    struct op {
-      static inline void initialise (T& accum) { accum = T(0); }
-      static inline void reduce (T& accum, const T& val, const CCTK_REAL weight) { if (weight>0) accum += weight*val*val; }
-      static inline void finalise (T& accum, const T& cnt) { }
-    };
-    MPI_Op mpi_op () const { return MPI_SUM; }
-  };
-  
-  struct average : reduction {
-    average () { }
-    ared thered () const { return do_average; }
+  struct mean : reduction {
+    ared thered () const { return do_mean; }
     bool uses_cnt () const { return true; }
     template<class T>
     struct op {
-      static inline void initialise (T& accum) { accum = T(0); }
-      static inline void reduce (T& accum, const T& val, const CCTK_REAL weight) { if (weight>0) accum += weight*val; }
+      static inline void initialise (T& accum) { accum = 0; }
+      static inline void reduce (T& accum, const T& val) { accum += val; }
       static inline void finalise (T& accum, const T& cnt) { accum /= cnt; }
     };
     MPI_Op mpi_op () const { return MPI_SUM; }
   };
   
   struct norm1 : reduction {
-    norm1 () { }
     ared thered () const { return do_norm1; }
     bool uses_cnt () const { return true; }
     template<class T>
     struct op {
-      static inline void initialise (T& accum) { accum = T(0); }
-      static inline void reduce (T& accum, const T& val, const CCTK_REAL weight) { if (weight>0) accum += weight*abs(val); }
+      static inline void initialise (T& accum) { accum = 0; }
+      static inline void reduce (T& accum, const T& val) { accum += abs(val); }
       static inline void finalise (T& accum, const T& cnt) { accum /= cnt; }
     };
     MPI_Op mpi_op () const { return MPI_SUM; }
   };
   
   struct norm2 : reduction {
-    norm2 () { }
     ared thered () const { return do_norm2; }
     bool uses_cnt () const { return true; }
     template<class T>
     struct op {
-      static inline void initialise (T& accum) { accum = T(0); }
-      static inline void reduce (T& accum, const T& val, const CCTK_REAL weight) { if (weight>0) accum += weight*abs(val)*abs(val); }
-      static inline void finalise (T& accum, const T& cnt) { accum = mysqrt(accum / cnt); }
+      static inline void initialise (T& accum) { accum = 0; }
+      static inline void reduce (T& accum, const T& val) { accum += abs(val)*abs(val); }
+      static inline void finalise (T& accum, const T& cnt) { accum = sqrt(accum / cnt); }
     };
     MPI_Op mpi_op () const { return MPI_SUM; }
   };
+  template<> inline void norm2::op<char>         ::finalise (char         & accum, const char         & cnt) { accum = sqrt((double)accum / cnt); }
+  template<> inline void norm2::op<signed char>  ::finalise (signed char  & accum, const signed char  & cnt) { accum = sqrt((double)accum / cnt); }
+  template<> inline void norm2::op<unsigned char>::finalise (unsigned char& accum, const unsigned char& cnt) { accum = sqrt((double)accum / cnt); }
+  template<> inline void norm2::op<short>        ::finalise (short        & accum, const short        & cnt) { accum = sqrt((double)accum / cnt); }
+  template<> inline void norm2::op<int>          ::finalise (int          & accum, const int          & cnt) { accum = sqrt((double)accum / cnt); }
+  template<> inline void norm2::op<long>         ::finalise (long         & accum, const long         & cnt) { accum = sqrt((double)accum / cnt); }
+  template<> inline void norm2::op<long long>    ::finalise (long long    & accum, const long long    & cnt) { accum = sqrt((double)accum / cnt); }
   
   struct norm_inf : reduction {
-    norm_inf () { }
     ared thered () const { return do_norm_inf; }
     bool uses_cnt () const { return false; }
     template<class T>
     struct op {
-      static inline void initialise (T& accum) { accum = T(0); }
-      static inline void reduce (T& accum, const T& val, const CCTK_REAL weight) { if (weight>0) accum = mymax(accum, T(abs(val))); }
+      static inline void initialise (T& accum) { accum = 0; }
+      static inline void reduce (T& accum, const T& val) { accum = max(accum, abs(val)); }
       static inline void finalise (T& accum, const T& cnt) { }
     };
-    MPI_Op mpi_op () const { return MPI_MAX; }
+    MPI_Op mpi_op () const { return MPI_SUM; }
   };
   
   
@@ -410,41 +193,22 @@ namespace CarpetReduce {
   void initialise (void* const outval, void* const cnt)
   {
     OP::initialise (*(T*)outval);
-    *(T*)cnt = T(0);
+    *(T*)cnt = 0;
   }
   
   template<class T,class OP>
-  void reduce (const int* const lsh, const int* const bbox,
-               const int* const nghostzones,
-	       vector<const void*> const& inarrays,
-               vector<CCTK_REAL> const& tfacs,
-               void* const outval, void* const cnt,
-               const CCTK_REAL* const weight, const CCTK_REAL levfac)
+  void reduce (const int* const dims, const int* const nghostzones,
+	       const void* const inarray, void* const outval, void* const cnt)
   {
-    for (size_t tl=0; tl<inarrays.size(); ++tl) {
-      assert (inarrays.at(tl));
-    }
-    assert (tfacs.size() == inarrays.size());
+    const T *myinarray = (const T*)inarray;
     T myoutval = *(T*)outval;
     T mycnt = *(T*)cnt;
-    vect<int,dim> imin, imax;
-    for (int d=0; d<dim; ++d) {
-      imin[d] = bbox[2*d  ] ? 0      : nghostzones[d];
-      imax[d] = bbox[2*d+1] ? lsh[d] : lsh[d] - nghostzones[d];
-    }
-    assert (dim==3);
-    for (int k=imin[2]; k<imax[2]; ++k) {
-      for (int j=imin[1]; j<imax[1]; ++j) {
-        for (int i=imin[0]; i<imax[0]; ++i) {
-	  const int index = i + lsh[0] * (j + lsh[1] * k);
-          CCTK_REAL const w = (weight ? weight[index] : 1.0) * levfac;
-          T myinval = T(0);
-          for (size_t tl=0; tl<inarrays.size(); ++tl) {
-            myinval += (static_cast<const T*>(inarrays.at(tl))[index]
-                        * tfacs.at(tl));
-          }
-	  OP::reduce (myoutval, myinval, w);
-	  mycnt += w;
+    for (int k=nghostzones[2]; k<dims[2]-nghostzones[2]; ++k) {
+      for (int j=nghostzones[1]; j<dims[1]-nghostzones[1]; ++j) {
+	for (int i=nghostzones[0]; i<dims[0]-nghostzones[0]; ++i) {
+	  const int index = i + dims[0] * (j + dims[1] * k);
+	  OP::reduce (myoutval, myinarray[index]);
+	  ++mycnt;
 	}
       }
     }
@@ -455,7 +219,7 @@ namespace CarpetReduce {
   template<class T,class OP>
   void finalise (void* const outval, const void* const cnt)
   {
-    OP::finalise (*(T*)outval, *(const T*)cnt);
+    OP::finalise (*(T*)outval, *(T*)cnt);
   }
   
   
@@ -483,29 +247,25 @@ namespace CarpetReduce {
     for (int n=0; n<num_outvals; ++n) {
       
       switch (outtype) {
-#define INITIALISE(OP,S)                                                \
-      case do_##OP: {                                                   \
-        typedef typeconv<S>::goodtype T;                                \
-	initialise<T,OP::op<T> > (&((char*)myoutvals)[vartypesize*n],   \
-				  &((char*)mycounts )[vartypesize*n]);  \
-	break;                                                          \
-      }
+#define INITIALISE(OP,T)						\
+      case do_##OP:							\
+	initialise<T,OP::op<T> > (&((char*)myoutvals)[vartypesize*n],	\
+				  &((char*)mycounts)[vartypesize*n]);	\
+	break;
 #define TYPECASE(N,T)				\
       case N: {					\
 	switch (red->thered()) {		\
 	  INITIALISE(count,T);			\
 	  INITIALISE(minimum,T);		\
 	  INITIALISE(maximum,T);		\
-	  INITIALISE(product,T);		\
 	  INITIALISE(sum,T);			\
 	  INITIALISE(sum_abs,T);		\
-	  INITIALISE(sum_squared,T);		\
-	  INITIALISE(average,T);		\
+	  INITIALISE(mean,T);			\
 	  INITIALISE(norm1,T);			\
 	  INITIALISE(norm2,T);			\
 	  INITIALISE(norm_inf,T);		\
 	default:				\
-	  assert (0);				\
+	  abort();				\
 	}					\
 	break;					\
       }
@@ -513,7 +273,7 @@ namespace CarpetReduce {
 #undef TYPECASE
 #undef INITIALISE
       default:
-	assert (0);
+	abort();
       }
       
     } // for n
@@ -521,72 +281,14 @@ namespace CarpetReduce {
   
   
   
-  void Copy (const cGH* const cgh, const int proc,
-             const int lsize,
-             const int num_inarrays,
-             const void* const* const inarrays, const int intype,
-             const int num_outvals,
-             void* const myoutvals, const int outtype,
-             void* const mycounts)
-  {
-    assert (cgh);
-    
-    assert (proc == -1 || (proc>=0 && proc<CCTK_nProcs(cgh)));
-    
-    assert (lsize >= 0);
-    assert (num_outvals>=0);
-    
-    assert (num_inarrays>=0);
-    assert (num_inarrays * lsize == num_outvals);
-    assert (inarrays);
-    for (int n=0; n<num_inarrays; ++n) {
-      assert (inarrays[n]);
-    }
-    
-    assert (myoutvals);
-    assert (mycounts);
-    
-    assert (outtype == intype);
-    
-    for (int m=0; m<num_inarrays; ++m) {
-      for (int n=0; n<lsize; ++n) {
-        
-        switch (outtype) {
-#define COPY(S)                                                         \
-        {                                                               \
-          typedef typeconv<S>::goodtype T;                              \
-          ((T*)myoutvals)[n+lsize*m] = ((const T*)inarrays[m])[n];      \
-          ((T*)mycounts )[n+lsize*m] = T(1);                            \
-        }
-#define TYPECASE(N,T)                           \
-        case N: {                               \
-          COPY(T);                              \
-          break;                                \
-        }
-#include "Carpet/Carpet/src/typecase"
-#undef TYPECASE
-#undef COPY
-        default:
-          assert (0);
-        }
-        
-      } // for
-    } // for
-  }
-  
-  
-  
   void Reduce (const cGH* const cgh, const int proc,
-	       const int* const mylsh, const int* const mybbox,
-               const int* const mynghostzones,
+	       const int* const mydims, const int* const mynghostzones,
 	       const int num_inarrays,
-	       vector<const void* const*> const& inarrays,
-               vector<CCTK_REAL> const& tfacs, const int intype,
+	       const void* const* const inarrays, const int intype,
 	       const int num_outvals,
 	       void* const myoutvals, const int outtype,
 	       void* const mycounts,
-	       const reduction* const red,
-               CCTK_REAL const * const weight, CCTK_REAL const levfac)
+	       const reduction* const red)
   {
     assert (cgh);
     
@@ -596,17 +298,14 @@ namespace CarpetReduce {
     
     assert (num_inarrays>=0);
     assert (num_inarrays == num_outvals);
-    for (size_t tl=0; tl<inarrays.size(); ++tl) {
-      assert (inarrays.at(tl));
-      for (int n=0; n<num_inarrays; ++n) {
-        assert (inarrays.at(tl)[n]);
-      }
+    assert (inarrays);
+    for (int n=0; n<num_inarrays; ++n) {
+      assert (inarrays[n]);
     }
-    assert (tfacs.size() == inarrays.size());
     
     for (int d=0; d<dim; ++d) {
-      assert (mylsh[d]>=0);
-      assert (mynghostzones[d]>=0 && 2*mynghostzones[d]<=mylsh[d]);
+      assert (mydims[d]>=0);
+      assert (mynghostzones[d]>=0 && 2*mynghostzones[d]<=mydims[d]);
     }
     
     const int vartypesize = CCTK_VarTypeSize(outtype);
@@ -615,43 +314,29 @@ namespace CarpetReduce {
     assert (myoutvals);
     assert (mycounts);
     
-    assert (outtype == intype);
-    
-    vector<const void*> myinarrays(inarrays.size());
-    
     for (int n=0; n<num_outvals; ++n) {
       
-      for (size_t tl=0; tl<inarrays.size(); ++tl) {
-        myinarrays.at(tl) = inarrays.at(tl)[n];
-      }
-      
       switch (outtype) {
-#define REDUCE(OP,S)                                                    \
-      case do_##OP: {                                                   \
-        typedef typeconv<S>::goodtype T;                                \
-        reduce<T,OP::op<T> > (mylsh, mybbox, mynghostzones,             \
-                              myinarrays, tfacs,                        \
-                              &((char*)myoutvals)[vartypesize*n],       \
-                              &((char*)mycounts )[vartypesize*n],       \
-                              weight, levfac);                          \
-        break;                                                          \
-      }
+#define REDUCE(OP,T)							\
+	case do_##OP:							\
+	  reduce<T,OP::op<T> > (mydims, mynghostzones, inarrays[n],	\
+				&((char*)myoutvals)[vartypesize*n],	\
+				&((char*)mycounts)[vartypesize*n]);	\
+	  break;
 #define TYPECASE(N,T)				\
       case N: {					\
 	switch (red->thered()) {		\
 	  REDUCE(count,T);			\
 	  REDUCE(minimum,T);			\
 	  REDUCE(maximum,T);			\
-	  REDUCE(product,T);			\
 	  REDUCE(sum,T);			\
 	  REDUCE(sum_abs,T);			\
-	  REDUCE(sum_squared,T);		\
-	  REDUCE(average,T);			\
+	  REDUCE(mean,T);			\
 	  REDUCE(norm1,T);			\
 	  REDUCE(norm2,T);			\
 	  REDUCE(norm_inf,T);			\
 	default:				\
-	  assert (0);				\
+	  abort();				\
 	}					\
 	break;					\
       }
@@ -659,7 +344,7 @@ namespace CarpetReduce {
 #undef TYPECASE
 #undef REDUCE
       default:
-	assert (0);
+	abort();
       }
       
     } // for n
@@ -687,29 +372,27 @@ namespace CarpetReduce {
     assert (myoutvals);
     assert (mycounts);
     
-    vector<char> counts;
+    char *counts = 0;
     if (proc==-1 || proc==CCTK_MyProc(cgh)) {
-      counts.resize(vartypesize * num_outvals);
+      counts = new char [vartypesize * num_outvals];
     }
     
-    const MPI_Datatype mpitype = CarpetSimpleMPIDatatype(outtype);
-    const int mpilength = CarpetSimpleMPIDatatypeLength(outtype);
     if (proc == -1) {
-      MPI_Allreduce ((void*)myoutvals, outvals, mpilength*num_outvals,
-		     mpitype, red->mpi_op(),
+      MPI_Allreduce ((void*)myoutvals, outvals, num_outvals,
+		     CarpetMPIDatatype(outtype), red->mpi_op(),
 		     CarpetMPIComm());
       if (red->uses_cnt()) {
-	MPI_Allreduce ((void*)mycounts, &counts[0], num_outvals*mpilength,
-		       mpitype, MPI_SUM,
+	MPI_Allreduce ((void*)mycounts, counts, num_outvals,
+		       CarpetMPIDatatype(outtype), MPI_SUM,
 		       CarpetMPIComm());
       }
     } else {
-      MPI_Reduce ((void*)myoutvals, outvals, num_outvals*mpilength,
-		  mpitype, red->mpi_op(),
+      MPI_Reduce ((void*)myoutvals, outvals, num_outvals,
+		  CarpetMPIDatatype(outtype), red->mpi_op(),
 		  proc, CarpetMPIComm());
       if (red->uses_cnt()) {
-	MPI_Reduce ((void*)mycounts, &counts[0], num_outvals*mpilength,
-		    mpitype, MPI_SUM,
+	MPI_Reduce ((void*)mycounts, counts, num_outvals,
+		    CarpetMPIDatatype(outtype), MPI_SUM,
 		    proc, CarpetMPIComm());
       }
     }
@@ -719,32 +402,28 @@ namespace CarpetReduce {
       for (int n=0; n<num_outvals; ++n) {
 	
 	assert (outvals);
-	assert ((int)counts.size() == vartypesize * num_outvals);
+	assert (counts);
 	
 	switch (outtype) {
-#define FINALISE(OP,S)                                                  \
-	case do_##OP: {                                                 \
-          typedef typeconv<S>::goodtype T;                              \
-	  finalise<T,OP::op<T> > (&((char*)outvals)[vartypesize*n],     \
-                                  &        counts  [vartypesize*n]);    \
-	  break;                                                        \
-        }
+#define FINALISE(OP,T)							\
+	case do_##OP:							\
+	  finalise<T,OP::op<T> >					\
+	    (&((char*)outvals)[vartypesize*n], &counts[vartypesize*n]);	\
+	  break;
 #define TYPECASE(N,T)				\
 	case N: {				\
 	  switch (red->thered()) {		\
 	    FINALISE(count,T);			\
 	    FINALISE(minimum,T);		\
 	    FINALISE(maximum,T);		\
-	    FINALISE(product,T);		\
 	    FINALISE(sum,T);			\
 	    FINALISE(sum_abs,T);		\
-	    FINALISE(sum_squared,T);		\
-	    FINALISE(average,T);		\
+	    FINALISE(mean,T);			\
 	    FINALISE(norm1,T);			\
 	    FINALISE(norm2,T);			\
 	    FINALISE(norm_inf,T);		\
 	  default:				\
-	    assert (0);				\
+	    abort();				\
 	  }					\
 	  break;				\
 	}
@@ -752,11 +431,13 @@ namespace CarpetReduce {
 #undef TYPECASE
 #undef FINALISE
 	default:
-	  assert (0);
+	  abort();
 	}
 	
       } // for n
       
+      delete [] counts;
+      
     } // if
   }
   
@@ -778,6 +459,7 @@ namespace CarpetReduce {
     assert (outvals || (proc!=-1 && proc!=CCTK_MyProc(cgh)));
     
     assert (num_inarrays>=0);
+    assert (num_inarrays == num_outvals);
     assert (inarrays);
     for (int n=0; n<num_inarrays; ++n) {
       assert (inarrays[n]);
@@ -788,58 +470,31 @@ namespace CarpetReduce {
       assert (dims[d]>=0);
     }
     
-    int lsize = 1;
+    int mydims[dim], mynghostzones[dim];
     for (int d=0; d<num_dims; ++d) {
-      lsize *= dims[d];
-    }
-    
-    const bool do_local_reduction = num_outvals == 1;
-    
-    if (! do_local_reduction) {
-      assert (num_outvals == lsize);
-    }
-    
-    vect<int,dim> mylsh, mynghostzones;
-    vect<vect<int,2>,dim> mybbox;
-    for (int d=0; d<num_dims; ++d) {
-      mylsh[d] = dims[d];
-      mybbox[d][0] = 0;
-      mybbox[d][1] = 0;
+      mydims[d] = dims[d];
       mynghostzones[d] = 0;
     }
     for (int d=num_dims; d<dim; ++d) {
-      mylsh[d] = 1;
-      mybbox[d][0] = 0;
-      mybbox[d][1] = 0;
+      mydims[d] = 1;
       mynghostzones[d] = 0;
     }
     
-    vector<const void* const*> myinarrays(1);
-    vector<CCTK_REAL> tfacs(1);
-    myinarrays.at(0) = inarrays;
-    tfacs.at(0) = 1.0;
-    
     const int vartypesize = CCTK_VarTypeSize(outtype);
     assert (vartypesize>=0);
     
-    vector<char> myoutvals (vartypesize * num_inarrays * num_outvals);
-    vector<char> mycounts  (vartypesize * num_inarrays * num_outvals);
+    char *myoutvals = new char [vartypesize * num_outvals];
+    char *mycounts  = new char [vartypesize * num_outvals];
     
-    Initialise (cgh, proc, num_inarrays * num_outvals, &myoutvals[0], outtype,
-                &mycounts[0], red);
-    if (do_local_reduction) {
-      Reduce   (cgh, proc, &mylsh[0], &mybbox[0][0], &mynghostzones[0],
-                num_inarrays, myinarrays, tfacs, intype,
-                num_inarrays * num_outvals, &myoutvals[0], outtype,
-                &mycounts[0], red,
-                NULL, 1.0);
-    } else {
-      Copy     (cgh, proc, lsize, num_inarrays, inarrays, intype,
-                num_inarrays * num_outvals, &myoutvals[0], outtype,
-                &mycounts[0]);
-    }
-    Finalise   (cgh, proc, num_inarrays * num_outvals, outvals, outtype,
-                &myoutvals[0], &mycounts[0], red);
+    Initialise (cgh, proc, num_outvals, myoutvals, outtype, mycounts, red);
+    Reduce     (cgh, proc, mydims, mynghostzones,
+		num_inarrays, inarrays, intype,
+		num_outvals, myoutvals, outtype, mycounts, red);
+    Finalise   (cgh, proc, num_outvals, outvals, outtype, myoutvals, mycounts,
+		red);
+    
+    delete [] myoutvals;
+    delete [] mycounts;
     
     return 0;
   }
@@ -853,289 +508,128 @@ namespace CarpetReduce {
   {
     int ierr;
     
+    // global mode
+    if (component != -1) {
+      CCTK_WARN (0, "It is not possible to use a grid variable reduction operator in local mode");
+    }
+    assert (component == -1);
+    
     assert (cgh);
     
     assert (proc == -1 || (proc>=0 && proc<CCTK_nProcs(cgh)));
     
     assert (num_outvals>=0);
-    assert (num_outvals==1);
     assert (outvals || (proc!=-1 && proc!=CCTK_MyProc(cgh)));
     
     assert (num_invars>=0);
+    assert (num_invars == num_outvals);
     assert (invars);
     for (int n=0; n<num_invars; ++n) {
       assert (invars[n]>=0 && invars[n]<CCTK_NumVars());
     }
     
-    if (num_invars==0) return 0;
+    if (num_outvals==0) return 0;
     
-    assert (num_invars>0);
     const int vi = invars[0];
     assert (vi>=0 && vi<CCTK_NumVars());
     
-    const int grpdim = CCTK_GroupDimFromVarI(vi);
-    assert (grpdim>=0 && grpdim<=dim);
-    for (int n=0; n<num_invars; ++n) {
-      assert (CCTK_GroupDimFromVarI(invars[n]) == grpdim);
-    }
-    
-    const int intype = CCTK_VarTypeI(vi);
+    const int num_dims = CCTK_GroupDimFromVarI(vi);
+    assert (num_dims>=0 && num_dims<=dim);
     for (int n=0; n<num_invars; ++n) {
-      assert (CCTK_VarTypeI(invars[n]) == intype);
+      assert (CCTK_GroupDimFromVarI(invars[n]) == num_dims);
     }
     
     const int vartypesize = CCTK_VarTypeSize(outtype);
     assert (vartypesize>=0);
     
+    char *myoutvals = new char [vartypesize * num_outvals];
+    char *mycounts  = new char [vartypesize * num_outvals];
     
+    Initialise (cgh, proc, num_outvals, myoutvals, outtype, mycounts, red);
     
-    // meta mode
-    if (is_meta_mode()) {
-      CCTK_WARN (0, "Grid variable reductions are not possible in meta mode");
-    }
-    
-    bool const reduce_arrays = CCTK_GroupTypeFromVarI(vi) != CCTK_GF;
-    bool const want_global_mode = is_global_mode() && ! reduce_arrays;
-    bool const want_level_mode = is_level_mode() && ! reduce_arrays;
-    
-    for (int n=0; n<num_invars; ++n) {
-      if ((CCTK_GroupTypeFromVarI(invars[n]) != CCTK_GF) != reduce_arrays) {
-        CCTK_WARN (0, "Cannot (yet) reduce grid functions and grid arrays/scalars at the same time");
+    BEGIN_COMPONENT_LOOP(cgh) {
+      
+      int dims[dim], nghostzones[dim];
+      ierr = CCTK_GrouplshVI(cgh, num_dims, dims, vi);
+      assert (!ierr);
+      ierr = CCTK_GroupnghostzonesVI(cgh, num_dims, nghostzones, vi);
+      assert (!ierr);
+      for (int d=0; d<num_dims; ++d) {
+	assert (dims[d]>=0);
+	assert (nghostzones[d]>=0 && 2*nghostzones[d]<=dims[d]);
       }
-    }
-    
-    // Ensure that all maps have the same number of refinement levels
-    for (int m=0; m<(int)vhh.size(); ++m) {
-      assert (vhh.at(m)->reflevels() == vhh.at(0)->reflevels());
-    }
-    int const minrl = reduce_arrays ? 0 : want_global_mode ? 0                      : reflevel;
-    int const maxrl = reduce_arrays ? 1 : want_global_mode ? vhh.at(0)->reflevels() : reflevel+1;
-    int const minm = reduce_arrays ? 0 : want_global_mode || want_level_mode ? 0    : Carpet::map;
-    int const maxm = reduce_arrays ? 1 : want_global_mode || want_level_mode ? maps : Carpet::map+1;
-    
-    
-    
-    // Find the time interpolation order
-    int partype;
-    void const * const parptr
-      = CCTK_ParameterGet ("prolongation_order_time", "Carpet", &partype);
-    assert (parptr);
-    assert (partype == PARAMETER_INTEGER);
-    int const prolongation_order_time = * (CCTK_INT const *) parptr;
-    
-    CCTK_REAL const current_time = cgh->cctk_time / cgh->cctk_delta_time;
-    
-    
-    
-    vector<char> myoutvals (vartypesize * num_invars * num_outvals);
-    vector<char> mycounts  (vartypesize * num_invars * num_outvals);
-    
-    Initialise (cgh, proc, num_invars * num_outvals, &myoutvals[0], outtype,
-                &mycounts[0], red);
+      
+      const void** const inarrays = new (const void*) [num_invars];
+      for (int n=0; n<num_invars; ++n) {
+	inarrays[n] = CCTK_VarDataPtrI(cgh, 0, invars[n]);
+	assert (inarrays[n]);
+      }
+      
+      const int intype = CCTK_VarTypeI(vi);
+      for (int n=0; n<num_invars; ++n) {
+	assert (CCTK_VarTypeI(invars[n]) == intype);
+      }
+      
+      int mydims[dim], mynghostzones[dim];
+      for (int d=0; d<num_dims; ++d) {
+	mydims[d] = dims[d];
+	mynghostzones[d] = nghostzones[d];
+      }
+      for (int d=num_dims; d<dim; ++d) {
+	mydims[d] = 1;
+	mynghostzones[d] = 0;
+      }
+      
+      Reduce (cgh, proc, mydims, mynghostzones, num_invars, inarrays, intype,
+	      num_outvals, myoutvals, outtype, mycounts, red);
+      
+      delete [] inarrays;
+      
+    } END_COMPONENT_LOOP(cgh);
     
-    BEGIN_GLOBAL_MODE(cgh) {
-      for (int rl=minrl; rl<maxrl; ++rl) {
-        enter_level_mode (const_cast<cGH*>(cgh), rl);
-        
-        
-        
-        // Number of necessary time levels
-        CCTK_REAL const level_time = cgh->cctk_time / cgh->cctk_delta_time;
-        bool need_time_interp
-          = (! reduce_arrays
-             && (fabs(current_time - level_time)
-                 > 1e-12 * (fabs(level_time) + fabs(current_time)
-                            + fabs(cgh->cctk_delta_time))));
-        assert (! (! want_global_mode && need_time_interp));
-        assert (! (reduce_arrays && need_time_interp));
-        int num_tl = need_time_interp ? prolongation_order_time + 1 : 1;
-        
-        // Are there enought time levels?
-        if (need_time_interp) {
-          
-          if (CCTK_ActiveTimeLevelsVI(cgh, vi) < num_tl) {
-            static vector<bool> have_warned;
-            if (have_warned.empty()) {
-              have_warned.resize (CCTK_NumVars(), false);
-            }
-            if (! have_warned.at(vi)) {
-              char * const fullname = CCTK_FullName(vi);
-              CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                          "Grid function \"%s\" has only %d active time levels on refinement level %d; this is not enough for time interpolation.  Using the current time level instead",
-                          fullname, CCTK_ActiveTimeLevelsVI(cgh, vi), reflevel);
-              free (fullname);
-              have_warned.at(vi) = true;
-            }
-            
-            // fall back
-            need_time_interp = false;
-            num_tl = 1;
-          }
-          
-        }
-        
-        vector<CCTK_REAL> tfacs(num_tl);
-        
-        // Interpolate in time, if necessary
-        if (need_time_interp) {
-          
-          // Get interpolation times
-          CCTK_REAL const time = current_time;
-          vector<CCTK_REAL> times(num_tl);
-          for (int tl=0; tl<num_tl; ++tl) {
-            times.at(tl) = vtt.at(0)->time (-tl, reflevel, mglevel);
-          }
-          
-          // Calculate interpolation weights
-          switch (num_tl) {
-          case 1:
-            // no interpolation
-            assert (fabs((time - times.at(0)) / fabs(time + times.at(0) + cgh->cctk_delta_time)) < 1e-12);
-            tfacs.at(0) = 1.0;
-            break;
-          case 2:
-            // linear (2-point) interpolation
-            tfacs.at(0) = (time - times.at(1)) / (times.at(0) - times.at(1));
-            tfacs.at(1) = (time - times.at(0)) / (times.at(1) - times.at(0));
-            break;
-          case 3:
-            // quadratic (3-point) interpolation
-            tfacs.at(0) = (time - times.at(1)) * (time - times.at(2)) / ((times.at(0) - times.at(1)) * (times.at(0) - times.at(2)));
-            tfacs.at(1) = (time - times.at(0)) * (time - times.at(2)) / ((times.at(1) - times.at(0)) * (times.at(1) - times.at(2)));
-            tfacs.at(2) = (time - times.at(0)) * (time - times.at(1)) / ((times.at(2) - times.at(0)) * (times.at(2) - times.at(1)));
-            break;
-          default:
-            assert (0);
-          }
-          
-        } else { // if ! need_time_interp
-          
-          assert (num_tl == 1);
-          tfacs.at(0) = 1;
-          
-        } // if ! need_time_interp
-        
-        
-        
-        for (int m=minm; m<maxm; ++m) {
-          enter_singlemap_mode (const_cast<cGH*>(cgh), m);
-          BEGIN_LOCAL_COMPONENT_LOOP(cgh, reduce_arrays ? CCTK_ARRAY : CCTK_GF) {
-            
-            
-            
-            assert (grpdim<=dim);
-            int lsh[dim], bbox[2*dim], nghostzones[dim];
-            ierr = CCTK_GrouplshVI(cgh, grpdim, lsh, vi);
-            assert (!ierr);
-            ierr = CCTK_GroupbboxVI(cgh, 2*grpdim, bbox, vi);
-            assert (!ierr);
-            ierr = CCTK_GroupnghostzonesVI(cgh, grpdim, nghostzones, vi);
-            assert (!ierr);
-            for (int d=0; d<grpdim; ++d) {
-              assert (lsh[d]>=0);
-              assert (nghostzones[d]>=0 && 2*nghostzones[d]<=lsh[d]);
-            }
-            
-            vect<int,dim> mylsh, mynghostzones;
-            vect<vect<int,2>,dim> mybbox;
-            for (int d=0; d<grpdim; ++d) {
-              mylsh[d] = lsh[d];
-              mybbox[d][0] = bbox[2*d  ];
-              mybbox[d][1] = bbox[2*d+1];
-              mynghostzones[d] = nghostzones[d];
-            }
-            for (int d=grpdim; d<dim; ++d) {
-              mylsh[d] = 1;
-              mybbox[d][0] = 0;
-              mybbox[d][1] = 0;
-              mynghostzones[d] = 0;
-            }
-            
-            
-            
-            CCTK_REAL const * weight;
-            CCTK_REAL levfac;
-            if (want_global_mode) {
-              weight = (static_cast<CCTK_REAL const *>
-                        (CCTK_VarDataPtr (cgh, 0, "CarpetReduce::weight")));
-              assert (weight);
-              levfac = pow(CCTK_REAL(reflevelfact), -grpdim);
-            } else {
-              weight = NULL;
-              levfac = 1.0;
-            }
-            
-            vector<vector<const void*> > myinarrays (num_tl);
-            vector<const void* const*> inarrays (num_tl);
-            for (int tl=0; tl<num_tl; ++tl) {
-              myinarrays.at(tl).resize (num_invars);
-              for (int n=0; n<num_invars; ++n) {
-                myinarrays.at(tl).at(n) = CCTK_VarDataPtrI(cgh, tl, invars[n]);
-                assert (myinarrays.at(tl).at(n));
-              }
-              inarrays.at(tl) = &myinarrays.at(tl).at(0);
-            }
-            
-            
-            
-            Reduce (cgh, proc, &mylsh[0], &mybbox[0][0], &mynghostzones[0],
-                    num_invars, inarrays, tfacs, intype,
-                    num_invars * num_outvals, &myoutvals[0], outtype,
-                    &mycounts[0], red,
-                    weight, levfac);
-            
-            
-            
-          } END_LOCAL_COMPONENT_LOOP;
-          leave_singlemap_mode (const_cast<cGH*>(cgh));
-        } // for m
-        
-        leave_level_mode (const_cast<cGH*>(cgh));
-      } // for rl
-    } END_GLOBAL_MODE;
+    Finalise (cgh, proc, num_outvals, outvals, outtype, myoutvals, mycounts,
+	      red);
     
-    Finalise (cgh, proc, num_invars * num_outvals, outvals, outtype,
-              &myoutvals[0], &mycounts[0], red);
+    delete [] myoutvals;
+    delete [] mycounts;
     
     return 0;
   }
   
   
   
-#define REDUCTION(OP)                                                     \
-  int OP##_arrays (const cGH * const cgh, const int proc,                 \
-		   const int num_dims, const int * const dims,            \
-		   const int num_inarrays,                                \
-		   const void * const * const inarrays, const int intype, \
-		   const int num_outvals,                                 \
-		   void * const outvals, const int outtype)               \
-  {                                                                       \
-    const OP red;                                                         \
-    return ReduceArrays                                                   \
-      (cgh, proc, num_dims, dims,                                         \
-       num_inarrays, inarrays, intype, num_outvals, outvals, outtype,     \
-       &red);                                                             \
-  }                                                                       \
-                                                                          \
-  int OP##_GVs (const cGH * const cgh, const int proc,                    \
-	        const int num_outvals,                                    \
-	        const int outtype, void * const outvals,                  \
-	        const int num_invars, const int * const invars)           \
-  {                                                                       \
-    const OP red;                                                         \
-    return ReduceGVs (cgh, proc,                                          \
-		      num_outvals, outtype, outvals, num_invars, invars,  \
-		      &red);                                              \
+#define REDUCTION(OP)							 \
+  int OP##_arrays (cGH * const cgh, const int proc,			 \
+		   const int num_dims, int * const dims,		 \
+		   const int num_inarrays,				 \
+		   void ** const inarrays, const int intype,		 \
+		   const int num_outvals,				 \
+		   void * const outvals, const int outtype)		 \
+  {									 \
+    const OP red;							 \
+    return ReduceArrays							 \
+      (cgh, proc, num_dims, dims,					 \
+       num_inarrays, inarrays, intype, num_outvals, outvals, outtype,	 \
+       &red);								 \
+  }									 \
+  									 \
+  int OP##_GVs (cGH * const cgh, const int proc,			 \
+	        const int num_outvals,					 \
+	        const int outtype, void * const outvals,		 \
+	        const int num_invars, int * const invars)		 \
+  {									 \
+    const OP red;							 \
+    return ReduceGVs (cgh, proc,					 \
+		      num_outvals, outtype, outvals, num_invars, invars, \
+		      &red);						 \
   }
   
   REDUCTION(count);
   REDUCTION(minimum);
   REDUCTION(maximum);
-  REDUCTION(product);
   REDUCTION(sum);
   REDUCTION(sum_abs);
-  REDUCTION(sum_squared);
-  REDUCTION(average);
+  REDUCTION(mean);
   REDUCTION(norm1);
   REDUCTION(norm2);
   REDUCTION(norm_inf);
@@ -1144,31 +638,29 @@ namespace CarpetReduce {
   
   
   
-  void CarpetReduceStartup ()
+  int CarpetReduceStartup ()
   {
-    CCTK_RegisterReductionOperator (count_GVs,       "count");
-    CCTK_RegisterReductionOperator (minimum_GVs,     "minimum");
-    CCTK_RegisterReductionOperator (maximum_GVs,     "maximum");
-    CCTK_RegisterReductionOperator (product_GVs,     "product");
-    CCTK_RegisterReductionOperator (sum_GVs,         "sum");
-    CCTK_RegisterReductionOperator (sum_abs_GVs,     "sum_abs");
-    CCTK_RegisterReductionOperator (sum_squared_GVs, "sum_squared");
-    CCTK_RegisterReductionOperator (average_GVs,     "average");
-    CCTK_RegisterReductionOperator (norm1_GVs,       "norm1");
-    CCTK_RegisterReductionOperator (norm2_GVs,       "norm2");
-    CCTK_RegisterReductionOperator (norm_inf_GVs,    "norm_inf");
+    CCTK_RegisterReductionOperator (count_GVs,    "count");
+    CCTK_RegisterReductionOperator (minimum_GVs,  "minimum");
+    CCTK_RegisterReductionOperator (maximum_GVs,  "maximum");
+    CCTK_RegisterReductionOperator (sum_GVs,      "sum");
+    CCTK_RegisterReductionOperator (sum_abs_GVs,  "sum_abs");
+    CCTK_RegisterReductionOperator (mean_GVs,     "mean");
+    CCTK_RegisterReductionOperator (norm1_GVs,    "norm1");
+    CCTK_RegisterReductionOperator (norm2_GVs,    "norm2");
+    CCTK_RegisterReductionOperator (norm_inf_GVs, "norm_inf");
+    
+    CCTK_RegisterReductionArrayOperator (count_arrays,    "count");
+    CCTK_RegisterReductionArrayOperator (minimum_arrays,  "minimum");
+    CCTK_RegisterReductionArrayOperator (maximum_arrays,  "maximum");
+    CCTK_RegisterReductionArrayOperator (sum_arrays,      "sum");
+    CCTK_RegisterReductionArrayOperator (sum_abs_arrays,  "sum_abs");
+    CCTK_RegisterReductionArrayOperator (mean_arrays,     "mean");
+    CCTK_RegisterReductionArrayOperator (norm1_arrays,    "norm1");
+    CCTK_RegisterReductionArrayOperator (norm2_arrays,    "norm2");
+    CCTK_RegisterReductionArrayOperator (norm_inf_arrays, "norm_inf");
     
-    CCTK_RegisterReductionArrayOperator (count_arrays,       "count");
-    CCTK_RegisterReductionArrayOperator (minimum_arrays,     "minimum");
-    CCTK_RegisterReductionArrayOperator (maximum_arrays,     "maximum");
-    CCTK_RegisterReductionArrayOperator (product_arrays,     "product");
-    CCTK_RegisterReductionArrayOperator (sum_arrays,         "sum");
-    CCTK_RegisterReductionArrayOperator (sum_abs_arrays,     "sum_abs");
-    CCTK_RegisterReductionArrayOperator (sum_squared_arrays, "sum_squared");
-    CCTK_RegisterReductionArrayOperator (average_arrays,     "average");
-    CCTK_RegisterReductionArrayOperator (norm1_arrays,       "norm1");
-    CCTK_RegisterReductionArrayOperator (norm2_arrays,       "norm2");
-    CCTK_RegisterReductionArrayOperator (norm_inf_arrays,    "norm_inf");
+    return 0;
   }
   
 } // namespace CarpetReduce
diff --git a/Carpet/CarpetReduce/src/reduce.h b/Carpet/CarpetReduce/src/reduce.h
index a30a5a735..9580fbfa1 100644
--- a/Carpet/CarpetReduce/src/reduce.h
+++ b/Carpet/CarpetReduce/src/reduce.h
@@ -1,19 +1,16 @@
-/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.h,v 1.3 2003/04/30 12:41:02 schnetter Exp $ */
+/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.h,v 1.1 2001/12/11 13:08:58 schnetter Exp $ */
 
 #ifndef CARPETREDUCE_H
 #define CARPETREDUCE_H
 
 #ifdef __cplusplus
-namespace CarpetReduce {
-  extern "C" {
+extern "C" {
 #endif
   
-    /* Scheduled functions */
-    void CarpetReduceStartup (void);
-    
+  int CarpetReduceStartup (void);
+  
 #ifdef __cplusplus
-  } /* extern "C" */
-} /* namespace CarpetReduce */
+} /* extern "C" */
 #endif
 
 #endif /* !defined(CARPETREDUCE_H) */
diff --git a/Carpet/CarpetReduce/src/reduce.hh b/Carpet/CarpetReduce/src/reduce.hh
index 856f2f823..b2508ee34 100644
--- a/Carpet/CarpetReduce/src/reduce.hh
+++ b/Carpet/CarpetReduce/src/reduce.hh
@@ -1,8 +1,12 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.hh,v 1.4 2002/09/01 14:52:28 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.hh,v 1.1 2001/12/11 13:08:59 schnetter Exp $
 
 #ifndef CARPETREDUCE_HH
 #define CARPETREDUCE_HH
 
+namespace CarpetReduce {
+  
 #include "reduce.h"
+  
+} // namespace CarpetReduce
 
 #endif // !defined(CARPETREDUCE_HH)