Parallelise loops via OpenMP.

Use C99 style variable declarations to simplify using OpenMP.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/InterpToArray/trunk@14 6ca9aeac-0e4f-0410-a746-fe4df63e9d0c

1 files changed, 287 insertions, 228 deletions

diff --git a/src/interp.c b/src/interp.c
index fa3ad16..d5d67b5 100644
--- a/src/interp.c
+++ b/src/interp.c
@@ -15,97 +15,80 @@ InterpToArray (CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
-  int interpolator;
-  int options_table;
-  int coord_handle;
   
-  CCTK_REAL * restrict coordsx;
-  CCTK_REAL * restrict coordsy;
-  CCTK_REAL * restrict coordsz;
-  CCTK_POINTER_TO_CONST coords[3];
-  CCTK_INT * restrict inputs;
-  CCTK_INT * restrict output_types;
-  CCTK_POINTER * restrict outputs;
-  CCTK_INT * restrict operation_codes;
-  CCTK_INT * restrict time_deriv_order;
   
-  int group;
-  cGroupDynamicData dyndata;
-  
-  int nvars;
-  int npoints;
-  
-  int n;
-  int i, j, k;
-  int ierr;
-  
-  
-  
-  interpolator = CCTK_InterpHandle (interpolator_name);
+  int const interpolator = CCTK_InterpHandle (interpolator_name);
   assert (interpolator >= 0);
   
-  options_table = Util_TableCreateFromString (interpolator_options);
+  int const options_table = Util_TableCreateFromString (interpolator_options);
   assert (options_table >= 0);
   
-  coord_handle = CCTK_CoordSystemHandle (interpolator_coordinates);
+  int const coord_handle = CCTK_CoordSystemHandle (interpolator_coordinates);
   assert (coord_handle >= 0);
   
   
   
   /* Scalars */
   {
-    nvars = nscalars;
+    int const nvars = nscalars;
     if (nvars > 0) {
-      npoints = 1;
+      int const npoints = 1;
       
-      coordsx = malloc (npoints * sizeof * coordsx);
+      CCTK_REAL * restrict const
+        coordsx = malloc (npoints * sizeof * coordsx);
       assert (npoints==0 || coordsx);
-      coordsy = malloc (npoints * sizeof * coordsy);
+      CCTK_REAL * restrict const
+        coordsy = malloc (npoints * sizeof * coordsy);
       assert (npoints==0 || coordsy);
-      coordsz = malloc (npoints * sizeof * coordsz);
+      CCTK_REAL * restrict const
+        coordsz = malloc (npoints * sizeof * coordsz);
       assert (npoints==0 || coordsz);
+      CCTK_POINTER_TO_CONST coords[3];
       coords[0] = coordsx;
       coords[1] = coordsy;
       coords[2] = coordsz;
       
-      n = 0;
-      assert (n <= npoints);
-      coordsx[n] = scalar_x0;
-      coordsy[n] = scalar_y0;
-      coordsz[n] = scalar_z0;
-      ++n;
-      assert (n == npoints);
+      {
+        int const n = 0;
+        coordsx[n] = scalar_x0;
+        coordsy[n] = scalar_y0;
+        coordsz[n] = scalar_z0;
+      }
       
-      inputs = malloc (nvars * sizeof * inputs);
+      CCTK_INT * restrict const
+        inputs = malloc (nvars * sizeof * inputs);
       assert (inputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         inputs[n] = CCTK_VarIndex (scalar_vars[n]);
         if (inputs[n] < 0) {
           inputs[n] = -1;
         }
       }
       
-      output_types = malloc (nvars * sizeof * output_types);
+      CCTK_INT * restrict const
+        output_types = malloc (nvars * sizeof * output_types);
       assert (output_types);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         output_types[n] = CCTK_VARIABLE_REAL;
       }
       
-      outputs = malloc (nvars * sizeof * outputs);
+      CCTK_POINTER * restrict const outputs = malloc (nvars * sizeof * outputs);
       assert (outputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         outputs[n] = &scalars[npoints * n];
       }
       
-      ierr = CCTK_InterpGridArrays
-        (cctkGH, 3, interpolator, options_table, coord_handle,
-         npoints, CCTK_VARIABLE_REAL, coords,
-         nvars, inputs,
-         nvars, output_types, outputs);
-      assert (! ierr);
+      {
+        int const ierr = CCTK_InterpGridArrays
+          (cctkGH, 3, interpolator, options_table, coord_handle,
+           npoints, CCTK_VARIABLE_REAL, coords,
+           nvars, inputs,
+           nvars, output_types, outputs);
+        assert (! ierr);
+      }
       
       free (coordsx);
       free (coordsy);
@@ -120,34 +103,37 @@ InterpToArray (CCTK_ARGUMENTS)
   
   /* 1D Arrays */
   {
-    nvars = narrays1d;
+    int const nvars = narrays1d;
     if (nvars > 0) {
-      npoints = array1d_npoints_i;
+      int const npoints = array1d_npoints_i;
       
-      coordsx = malloc (npoints * sizeof * coordsx);
+      CCTK_REAL * restrict const
+        coordsx = malloc (npoints * sizeof * coordsx);
       assert (npoints==0 || coordsx);
-      coordsy = malloc (npoints * sizeof * coordsy);
+      CCTK_REAL * restrict const
+        coordsy = malloc (npoints * sizeof * coordsy);
       assert (npoints==0 || coordsy);
-      coordsz = malloc (npoints * sizeof * coordsz);
+      CCTK_REAL * restrict const
+        coordsz = malloc (npoints * sizeof * coordsz);
       assert (npoints==0 || coordsz);
+      CCTK_POINTER_TO_CONST coords[3];
       coords[0] = coordsx;
       coords[1] = coordsy;
       coords[2] = coordsz;
       
-      n = 0;
-      for (i=0; i<array1d_npoints_i; ++i) {
-        assert (n <= npoints);
+#pragma omp parallel for
+      for (int i=0; i<array1d_npoints_i; ++i) {
+        int const n = i;
         coordsx[n] = array1d_x0 + i * array1d_dx_i;
         coordsy[n] = array1d_y0 + i * array1d_dy_i;
         coordsz[n] = array1d_z0 + i * array1d_dz_i;
-        ++n;
       }
-      assert (n == npoints);
       
-      inputs = malloc (nvars * sizeof * inputs);
+      CCTK_INT * restrict const
+        inputs = malloc (nvars * sizeof * inputs);
       assert (inputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         inputs[n] = CCTK_VarIndex (array1d_vars[n]);
         assert (inputs[n] >= 0);
         if (inputs[n] < 0) {
@@ -155,54 +141,69 @@ InterpToArray (CCTK_ARGUMENTS)
         }
       }
       
-      output_types = malloc (nvars * sizeof * output_types);
+      CCTK_INT * restrict const
+        output_types = malloc (nvars * sizeof * output_types);
       assert (output_types);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         output_types[n] = CCTK_VARIABLE_REAL;
       }
       
-      outputs = malloc (nvars * sizeof * outputs);
+      CCTK_POINTER * restrict const
+        outputs = malloc (nvars * sizeof * outputs);
       assert (outputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         outputs[n] = &arrays1d[npoints * n];
       }
       
-      operation_codes = malloc (nvars * sizeof * operation_codes);
+      CCTK_INT * restrict const
+        operation_codes = malloc (nvars * sizeof * operation_codes);
       assert (operation_codes);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         operation_codes[n] = array1d_spacederivs[n];
       }
       
-      ierr = Util_TableSetIntArray
-        (options_table, nvars, operation_codes, "operation_codes");
-      assert (! ierr);
+      {
+        int const ierr = Util_TableSetIntArray
+          (options_table, nvars, operation_codes, "operation_codes");
+        assert (! ierr);
+      }
       
-      time_deriv_order = malloc (nvars * sizeof * time_deriv_order);
+      CCTK_INT * restrict const
+        time_deriv_order = malloc (nvars * sizeof * time_deriv_order);
       assert (time_deriv_order);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         time_deriv_order[n] = array1d_timederivs[n];
       }
       
-      ierr = Util_TableSetIntArray
-        (options_table, nvars, time_deriv_order, "time_deriv_order");
-      assert (! ierr);
+      {
+        int const ierr = Util_TableSetIntArray
+          (options_table, nvars, time_deriv_order, "time_deriv_order");
+        assert (! ierr);
+      }
       
-      ierr = CCTK_InterpGridArrays
-        (cctkGH, 3, interpolator, options_table, coord_handle,
-         npoints, CCTK_VARIABLE_REAL, coords,
-         nvars, inputs,
-         nvars, output_types, outputs);
-      assert (! ierr);
+      {
+        int const ierr = CCTK_InterpGridArrays
+          (cctkGH, 3, interpolator, options_table, coord_handle,
+           npoints, CCTK_VARIABLE_REAL, coords,
+           nvars, inputs,
+           nvars, output_types, outputs);
+        assert (! ierr);
+      }
       
-      ierr = Util_TableDeleteKey (options_table, "operation_codes");
-      assert (! ierr);
+      {
+        int const ierr = Util_TableDeleteKey (options_table, "operation_codes");
+        assert (! ierr);
+      }
       
-      ierr = Util_TableDeleteKey (options_table, "time_deriv_order");
-      assert (! ierr);
+      {
+        int const ierr =
+          Util_TableDeleteKey (options_table, "time_deriv_order");
+        assert (! ierr);
+      }
       
       free (coordsx);
       free (coordsy);
@@ -219,62 +220,69 @@ InterpToArray (CCTK_ARGUMENTS)
   
   /* 2D Arrays */
   {
-    nvars = narrays2d;
+    int const nvars = narrays2d;
     if (nvars > 0) {
-      npoints = array2d_npoints_i * array2d_npoints_j;
+      int const npoints = array2d_npoints_i * array2d_npoints_j;
       
-      coordsx = malloc (npoints * sizeof * coordsx);
+      CCTK_REAL * restrict const
+        coordsx = malloc (npoints * sizeof * coordsx);
       assert (npoints==0 || coordsx);
-      coordsy = malloc (npoints * sizeof * coordsy);
+      CCTK_REAL * restrict const
+        coordsy = malloc (npoints * sizeof * coordsy);
       assert (npoints==0 || coordsy);
-      coordsz = malloc (npoints * sizeof * coordsz);
+      CCTK_REAL * restrict const
+        coordsz = malloc (npoints * sizeof * coordsz);
       assert (npoints==0 || coordsz);
+      CCTK_POINTER_TO_CONST coords[3];
       coords[0] = coordsx;
       coords[1] = coordsy;
       coords[2] = coordsz;
       
-      n = 0;
-      for (j=0; j<array2d_npoints_j; ++j) {
-        for (i=0; i<array2d_npoints_i; ++i) {
-          assert (n <= npoints);
+#pragma omp parallel for
+      for (int j=0; j<array2d_npoints_j; ++j) {
+        for (int i=0; i<array2d_npoints_i; ++i) {
+          const int n = i + array2d_npoints_i * j;
           coordsx[n] = array2d_x0 + i * array2d_dx_i + j * array2d_dx_j;
           coordsy[n] = array2d_y0 + i * array2d_dy_i + j * array2d_dy_j;
           coordsz[n] = array2d_z0 + i * array2d_dz_i + j * array2d_dz_j;
-          ++n;
         }
       }
-      assert (n == npoints);
       
-      inputs = malloc (nvars * sizeof * inputs);
+      CCTK_INT * restrict const
+        inputs = malloc (nvars * sizeof * inputs);
       assert (inputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         inputs[n] = CCTK_VarIndex (array2d_vars[n]);
         if (inputs[n] < 0) {
           inputs[n] = -1;
         }
       }
       
-      output_types = malloc (nvars * sizeof * output_types);
+      CCTK_INT * restrict const
+        output_types = malloc (nvars * sizeof * output_types);
       assert (output_types);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         output_types[n] = CCTK_VARIABLE_REAL;
       }
       
-      outputs = malloc (nvars * sizeof * outputs);
+      CCTK_POINTER * restrict const
+        outputs = malloc (nvars * sizeof * outputs);
       assert (outputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         outputs[n] = &arrays2d[npoints * n];
       }
       
-      ierr = CCTK_InterpGridArrays
-        (cctkGH, 3, interpolator, options_table, coord_handle,
-         npoints, CCTK_VARIABLE_REAL, coords,
-         nvars, inputs,
-         nvars, output_types, outputs);
-      assert (! ierr);
+      {
+        int const ierr = CCTK_InterpGridArrays
+          (cctkGH, 3, interpolator, options_table, coord_handle,
+           npoints, CCTK_VARIABLE_REAL, coords,
+           nvars, inputs,
+           nvars, output_types, outputs);
+        assert (! ierr);
+      }
       
       free (coordsx);
       free (coordsy);
@@ -289,64 +297,72 @@ InterpToArray (CCTK_ARGUMENTS)
   
   /* 3D Arrays */
   {
-    nvars = narrays3d;
+    int const nvars = narrays3d;
     if (nvars > 0) {
-      npoints = array3d_npoints_i * array3d_npoints_j * array3d_npoints_k;
+      int const npoints =
+        array3d_npoints_i * array3d_npoints_j * array3d_npoints_k;
       
-      coordsx = malloc (npoints * sizeof * coordsx);
+      CCTK_REAL * restrict const
+        coordsx = malloc (npoints * sizeof * coordsx);
       assert (npoints==0 || coordsx);
-      coordsy = malloc (npoints * sizeof * coordsy);
+      CCTK_REAL * restrict const
+        coordsy = malloc (npoints * sizeof * coordsy);
       assert (npoints==0 || coordsy);
-      coordsz = malloc (npoints * sizeof * coordsz);
+      CCTK_REAL * restrict const
+        coordsz = malloc (npoints * sizeof * coordsz);
       assert (npoints==0 || coordsz);
+      CCTK_POINTER_TO_CONST coords[3];
       coords[0] = coordsx;
       coords[1] = coordsy;
       coords[2] = coordsz;
       
-      n = 0;
-      for (k=0; k<array3d_npoints_k; ++k) {
-        for (j=0; j<array3d_npoints_j; ++j) {
-          for (i=0; i<array3d_npoints_i; ++i) {
-            assert (n <= npoints);
+#pragma omp parallel for
+      for (int k=0; k<array3d_npoints_k; ++k) {
+        for (int j=0; j<array3d_npoints_j; ++j) {
+          for (int i=0; i<array3d_npoints_i; ++i) {
+            int const n = i + array3d_npoints_i * (j + array3d_npoints_j * k);
             coordsx[n] = array3d_x0 + i * array3d_dx_i + j * array3d_dx_j + k * array3d_dx_k;
             coordsy[n] = array3d_y0 + i * array3d_dy_i + j * array3d_dy_j + k * array3d_dy_k;
             coordsz[n] = array3d_z0 + i * array3d_dz_i + j * array3d_dz_j + k * array3d_dz_k;
-            ++n;
           }
         }
       }
-      assert (n == npoints);
       
-      inputs = malloc (nvars * sizeof * inputs);
+      CCTK_INT * restrict const
+        inputs = malloc (nvars * sizeof * inputs);
       assert (inputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         inputs[n] = CCTK_VarIndex (array3d_vars[n]);
         if (inputs[n] < 0) {
           inputs[n] = -1;
         }
       }
       
-      output_types = malloc (nvars * sizeof * output_types);
+      CCTK_INT * restrict const
+        output_types = malloc (nvars * sizeof * output_types);
       assert (output_types);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         output_types[n] = CCTK_VARIABLE_REAL;
       }
       
-      outputs = malloc (nvars * sizeof * outputs);
+      CCTK_POINTER * restrict const
+        outputs = malloc (nvars * sizeof * outputs);
       assert (outputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         outputs[n] = &arrays3d[npoints * n];
       }
       
-      ierr = CCTK_InterpGridArrays
-        (cctkGH, 3, interpolator, options_table, coord_handle,
-         npoints, CCTK_VARIABLE_REAL, coords,
-         nvars, inputs,
-         nvars, output_types, outputs);
-      assert (! ierr);
+      {
+        int const ierr = CCTK_InterpGridArrays
+          (cctkGH, 3, interpolator, options_table, coord_handle,
+           npoints, CCTK_VARIABLE_REAL, coords,
+           nvars, inputs,
+           nvars, output_types, outputs);
+        assert (! ierr);
+      }
       
       free (coordsx);
       free (coordsy);
@@ -361,42 +377,48 @@ InterpToArray (CCTK_ARGUMENTS)
   
   /* Parallel 1D Arrays */
   {
-    nvars = nparrays1d;
+    int const nvars = nparrays1d;
     if (nvars > 0) {
-      group = CCTK_GroupIndex ("InterpToArray::parrays1d");
+      int const group = CCTK_GroupIndex ("InterpToArray::parrays1d");
       assert (group >= 0);
-      ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
-      assert (! ierr);
+      cGroupDynamicData dyndata;
+      {
+        int const ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
+        assert (! ierr);
+      }
       assert (dyndata.dim == 1);
-      npoints = dyndata.lsh[0];
+      int const npoints = dyndata.lsh[0];
       
-      coordsx = malloc (npoints * sizeof * coordsx);
+      CCTK_REAL * restrict const
+        coordsx = malloc (npoints * sizeof * coordsx);
       assert (npoints==0 || coordsx);
-      coordsy = malloc (npoints * sizeof * coordsy);
+      CCTK_REAL * restrict const
+        coordsy = malloc (npoints * sizeof * coordsy);
       assert (npoints==0 || coordsy);
-      coordsz = malloc (npoints * sizeof * coordsz);
+      CCTK_REAL * restrict const
+        coordsz = malloc (npoints * sizeof * coordsz);
       assert (npoints==0 || coordsz);
+      CCTK_POINTER_TO_CONST coords[3];
       coords[0] = coordsx;
       coords[1] = coordsy;
       coords[2] = coordsz;
       
-      n = 0;
-      for (i=0; i<dyndata.lsh[0]; ++i) {
-        assert (n <= npoints);
+#pragma omp parallel for
+      for (int i=0; i<dyndata.lsh[0]; ++i) {
+        int const n = i;
         coordsx[n] = (parray1d_x0 +
                       (dyndata.lbnd[0] + i) * parray1d_dx_i);
         coordsy[n] = (parray1d_y0 +
                       (dyndata.lbnd[0] + i) * parray1d_dy_i);
         coordsz[n] = (parray1d_z0 +
                       (dyndata.lbnd[0] + i) * parray1d_dz_i);
-        ++n;
       }
-      assert (n == npoints);
       
-      inputs = malloc (nvars * sizeof * inputs);
+      CCTK_INT * restrict const
+        inputs = malloc (nvars * sizeof * inputs);
       assert (inputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         inputs[n] = CCTK_VarIndex (parray1d_vars[n]);
         assert (inputs[n] >= 0);
         if (inputs[n] < 0) {
@@ -404,54 +426,69 @@ InterpToArray (CCTK_ARGUMENTS)
         }
       }
       
-      output_types = malloc (nvars * sizeof * output_types);
+      CCTK_INT * restrict const
+        output_types = malloc (nvars * sizeof * output_types);
       assert (output_types);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         output_types[n] = CCTK_VARIABLE_REAL;
       }
       
-      outputs = malloc (nvars * sizeof * outputs);
+      CCTK_POINTER * restrict const
+        outputs = malloc (nvars * sizeof * outputs);
       assert (outputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         outputs[n] = &parrays1d[npoints * n];
       }
       
-      operation_codes = malloc (nvars * sizeof * operation_codes);
+      CCTK_INT * restrict const
+        operation_codes = malloc (nvars * sizeof * operation_codes);
       assert (operation_codes);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         operation_codes[n] = parray1d_spacederivs[n];
       }
       
-      ierr = Util_TableSetIntArray
-        (options_table, nvars, operation_codes, "operation_codes");
-      assert (! ierr);
+      {
+        int const ierr = Util_TableSetIntArray
+          (options_table, nvars, operation_codes, "operation_codes");
+        assert (! ierr);
+      }
       
-      time_deriv_order = malloc (nvars * sizeof * time_deriv_order);
+      CCTK_INT * restrict const
+        time_deriv_order = malloc (nvars * sizeof * time_deriv_order);
       assert (time_deriv_order);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         time_deriv_order[n] = parray1d_timederivs[n];
       }
       
-      ierr = Util_TableSetIntArray
-        (options_table, nvars, time_deriv_order, "time_deriv_order");
-      assert (! ierr);
+      {
+        int const ierr = Util_TableSetIntArray
+          (options_table, nvars, time_deriv_order, "time_deriv_order");
+        assert (! ierr);
+      }
       
-      ierr = CCTK_InterpGridArrays
-        (cctkGH, 3, interpolator, options_table, coord_handle,
-         npoints, CCTK_VARIABLE_REAL, coords,
-         nvars, inputs,
-         nvars, output_types, outputs);
-      assert (! ierr);
+      {
+        int const ierr = CCTK_InterpGridArrays
+          (cctkGH, 3, interpolator, options_table, coord_handle,
+           npoints, CCTK_VARIABLE_REAL, coords,
+           nvars, inputs,
+           nvars, output_types, outputs);
+        assert (! ierr);
+      }
       
-      ierr = Util_TableDeleteKey (options_table, "operation_codes");
-      assert (! ierr);
+      {
+        int const ierr = Util_TableDeleteKey (options_table, "operation_codes");
+        assert (! ierr);
+      }
       
-      ierr = Util_TableDeleteKey (options_table, "time_deriv_order");
-      assert (! ierr);
+      {
+        int const ierr =
+          Util_TableDeleteKey (options_table, "time_deriv_order");
+        assert (! ierr);
+      }
       
       free (coordsx);
       free (coordsy);
@@ -468,29 +505,36 @@ InterpToArray (CCTK_ARGUMENTS)
   
   /* Parallel 2D Arrays */
   {
-    nvars = nparrays2d;
+    int const nvars = nparrays2d;
     if (nvars > 0) {
-      group = CCTK_GroupIndex ("InterpToArray::parrays2d");
+      int const group = CCTK_GroupIndex ("InterpToArray::parrays2d");
       assert (group >= 0);
-      ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
-      assert (! ierr);
+      cGroupDynamicData dyndata;
+      {
+        int const ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
+        assert (! ierr);
+      }
       assert (dyndata.dim == 2);
-      npoints = dyndata.lsh[0] * dyndata.lsh[1];
+      int const npoints = dyndata.lsh[0] * dyndata.lsh[1];
       
-      coordsx = malloc (npoints * sizeof * coordsx);
+      CCTK_REAL * restrict const
+        coordsx = malloc (npoints * sizeof * coordsx);
       assert (npoints==0 || coordsx);
-      coordsy = malloc (npoints * sizeof * coordsy);
+      CCTK_REAL * restrict const
+        coordsy = malloc (npoints * sizeof * coordsy);
       assert (npoints==0 || coordsy);
-      coordsz = malloc (npoints * sizeof * coordsz);
+      CCTK_REAL * restrict const
+        coordsz = malloc (npoints * sizeof * coordsz);
       assert (npoints==0 || coordsz);
+      CCTK_POINTER_TO_CONST coords[3];
       coords[0] = coordsx;
       coords[1] = coordsy;
       coords[2] = coordsz;
       
-      n = 0;
-      for (j=0; j<dyndata.lsh[1]; ++j) {
-        for (i=0; i<dyndata.lsh[0]; ++i) {
-          assert (n <= npoints);
+#pragma omp parallel for
+      for (int j=0; j<dyndata.lsh[1]; ++j) {
+        for (int i=0; i<dyndata.lsh[0]; ++i) {
+          int const n = i + dyndata.lsh[0] * j;
           coordsx[n] = (parray2d_x0 +
                         (dyndata.lbnd[0] + i) * parray2d_dx_i +
                         (dyndata.lbnd[1] + j) * parray2d_dx_j);
@@ -500,41 +544,44 @@ InterpToArray (CCTK_ARGUMENTS)
           coordsz[n] = (parray2d_z0 +
                         (dyndata.lbnd[0] + i) * parray2d_dz_i +
                         (dyndata.lbnd[1] + j) * parray2d_dz_j);
-          ++n;
         }
       }
-      assert (n == npoints);
       
-      inputs = malloc (nvars * sizeof * inputs);
+      CCTK_INT * restrict const
+        inputs = malloc (nvars * sizeof * inputs);
       assert (inputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         inputs[n] = CCTK_VarIndex (parray2d_vars[n]);
         if (inputs[n] < 0) {
           inputs[n] = -1;
         }
       }
       
-      output_types = malloc (nvars * sizeof * output_types);
+      CCTK_INT * restrict const
+        output_types = malloc (nvars * sizeof * output_types);
       assert (output_types);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         output_types[n] = CCTK_VARIABLE_REAL;
       }
       
-      outputs = malloc (nvars * sizeof * outputs);
+      CCTK_POINTER * restrict const
+        outputs = malloc (nvars * sizeof * outputs);
       assert (outputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         outputs[n] = &parrays2d[npoints * n];
       }
       
-      ierr = CCTK_InterpGridArrays
-        (cctkGH, 3, interpolator, options_table, coord_handle,
-         npoints, CCTK_VARIABLE_REAL, coords,
-         nvars, inputs,
-         nvars, output_types, outputs);
-      assert (! ierr);
+      {
+        int const ierr = CCTK_InterpGridArrays
+          (cctkGH, 3, interpolator, options_table, coord_handle,
+           npoints, CCTK_VARIABLE_REAL, coords,
+           nvars, inputs,
+           nvars, output_types, outputs);
+        assert (! ierr);
+      }
       
       free (coordsx);
       free (coordsy);
@@ -549,30 +596,37 @@ InterpToArray (CCTK_ARGUMENTS)
   
   /* Parallel 3D Arrays */
   {
-    nvars = nparrays3d;
+    int const nvars = nparrays3d;
     if (nvars > 0) {
-      group = CCTK_GroupIndex ("InterpToArray::parrays3d");
+      int const group = CCTK_GroupIndex ("InterpToArray::parrays3d");
       assert (group >= 0);
-      ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
-      assert (! ierr);
+      cGroupDynamicData dyndata;
+      {
+        int const ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
+        assert (! ierr);
+      }
       assert (dyndata.dim == 3);
-      npoints = dyndata.lsh[0] * dyndata.lsh[1] * dyndata.lsh[2];
+      int const npoints = dyndata.lsh[0] * dyndata.lsh[1] * dyndata.lsh[2];
       
-      coordsx = malloc (npoints * sizeof * coordsx);
+      CCTK_REAL * restrict const
+        coordsx = malloc (npoints * sizeof * coordsx);
       assert (npoints==0 || coordsx);
-      coordsy = malloc (npoints * sizeof * coordsy);
+      CCTK_REAL * restrict const
+        coordsy = malloc (npoints * sizeof * coordsy);
       assert (npoints==0 || coordsy);
-      coordsz = malloc (npoints * sizeof * coordsz);
+      CCTK_REAL * restrict const
+        coordsz = malloc (npoints * sizeof * coordsz);
       assert (npoints==0 || coordsz);
+      CCTK_POINTER_TO_CONST coords[3];
       coords[0] = coordsx;
       coords[1] = coordsy;
       coords[2] = coordsz;
       
-      n = 0;
-      for (k=0; k<dyndata.lsh[2]; ++k) {
-        for (j=0; j<dyndata.lsh[1]; ++j) {
-          for (i=0; i<dyndata.lsh[0]; ++i) {
-            assert (n <= npoints);
+#pragma omp parallel for
+      for (int k=0; k<dyndata.lsh[2]; ++k) {
+        for (int j=0; j<dyndata.lsh[1]; ++j) {
+          for (int i=0; i<dyndata.lsh[0]; ++i) {
+            int const n = i + dyndata.lsh[0] * (j + dyndata.lsh[1] * k);
             coordsx[n] = (parray3d_x0 +
                           (dyndata.lbnd[0] + i) * parray3d_dx_i +
                           (dyndata.lbnd[1] + j) * parray3d_dx_j +
@@ -585,42 +639,45 @@ InterpToArray (CCTK_ARGUMENTS)
                           (dyndata.lbnd[0] + i) * parray3d_dz_i +
                           (dyndata.lbnd[1] + j) * parray3d_dz_j +
                           (dyndata.lbnd[2] + k) * parray3d_dz_k);
-            ++n;
           }
         }
       }
-      assert (n == npoints);
       
-      inputs = malloc (nvars * sizeof * inputs);
+      CCTK_INT * restrict const
+        inputs = malloc (nvars * sizeof * inputs);
       assert (inputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         inputs[n] = CCTK_VarIndex (parray3d_vars[n]);
         if (inputs[n] < 0) {
           inputs[n] = -1;
         }
       }
       
-      output_types = malloc (nvars * sizeof * output_types);
+      CCTK_INT * restrict const
+        output_types = malloc (nvars * sizeof * output_types);
       assert (output_types);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         output_types[n] = CCTK_VARIABLE_REAL;
       }
       
-      outputs = malloc (nvars * sizeof * outputs);
+      CCTK_POINTER * restrict const
+        outputs = malloc (nvars * sizeof * outputs);
       assert (outputs);
       
-      for (n=0; n<nvars; ++n) {
+      for (int n=0; n<nvars; ++n) {
         outputs[n] = &parrays3d[npoints * n];
       }
       
-      ierr = CCTK_InterpGridArrays
-        (cctkGH, 3, interpolator, options_table, coord_handle,
-         npoints, CCTK_VARIABLE_REAL, coords,
-         nvars, inputs,
-         nvars, output_types, outputs);
-      assert (! ierr);
+      {
+        int const ierr = CCTK_InterpGridArrays
+          (cctkGH, 3, interpolator, options_table, coord_handle,
+           npoints, CCTK_VARIABLE_REAL, coords,
+           nvars, inputs,
+           nvars, output_types, outputs);
+        assert (! ierr);
+      }
       
       free (coordsx);
       free (coordsy);
@@ -633,6 +690,8 @@ InterpToArray (CCTK_ARGUMENTS)
   
   
   
-  ierr = Util_TableDestroy (options_table);
-  assert (! ierr);
+  {
+    int const ierr = Util_TableDestroy (options_table);
+    assert (! ierr);
+  }
 }