For the multiprocessor case, set the local interpolator status code to the

minimum of status codes from all the interpolation points a processor had
requested - provided that the local interpolator supports the feature of
returning such per-point status information (AEILocalInterp will do so soon).


git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHInterp/trunk@65 1c20744c-e24a-42ec-9533-f5004cb800e5

1 files changed, 103 insertions, 16 deletions

diff --git a/src/InterpGridArrays.c b/src/InterpGridArrays.c
index 23db30b..23089a4 100644
--- a/src/InterpGridArrays.c
+++ b/src/InterpGridArrays.c
@@ -197,7 +197,9 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
   pughInterpGH *myGH;
   CCTK_REAL **interp_coords_local;
   int N_points_local, N_types;
-  int j, point, proc, array, type, offset;
+  int j, point, proc, array, offset, type, table_handle;
+  CCTK_INT have_per_point_status, error_point_status;
+  CCTK_INT *per_point_status, *per_proc_status, *local_interp_status;
   char *msg;
   char **output_arrays_local;
   type_desc_t *this, *type_desc;
@@ -386,6 +388,11 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
     free (input_arrays);
     free (input_array_dims);
 
+    /* store this processor's interpolator status in the parameter table */
+    if (param_table_handle >= 0)
+    {
+      Util_TableSetInt(param_table_handle, retval, "local_interpolator_status");
+    }
     return (retval);
   }
 
@@ -828,6 +835,19 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
     free (coords);
   }
 
+  /* prepare the per_proc_status[_p_] status array which will contain the
+     per-point status codes from CCTK_InterpLocalUniform() reduced over all
+     points requested by processor _p_
+     Since the per_point_status is an optional return argument which may not
+     be supported by the local interpolator, we allocate the per_proc_status[]
+     array to be all zeros as a default status. */
+  have_per_point_status = 0;
+  per_proc_status = malloc ((2*pughGH->nprocs + N_points_local)
+                            * sizeof (CCTK_INT));
+  memset (per_proc_status, 0, pughGH->nprocs * sizeof (CCTK_INT));
+  local_interp_status = per_proc_status + 1*pughGH->nprocs;
+  per_point_status    = per_proc_status + 2*pughGH->nprocs;
+
   /* allocate intermediate output arrays for local interpolation */
   output_arrays_local = calloc (N_output_arrays, sizeof (void *));
   if (N_points_local > 0)
@@ -837,27 +857,89 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
       this = type_desc + output_array_types[array];
       output_arrays_local[array] = malloc (N_points_local * this->vtypesize);
     }
+  }
+  else
+  {
+    /* this indicates a query of the local interpolator whether it will
+       return per-point status information */
+    per_point_status = NULL;
+  }
 
-    /* now call the local interpolator for all local points and store
-       the results in the intermediate local output arrays */
-    retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle,
-                                      param_table_handle,
-                                      origin_local, delta, N_points_local,
-                                      interp_coords_type,
-                                      (const void **) interp_coords_local,
-                                      N_input_arrays, input_array_dims,
-                                      input_array_types, input_arrays,
-                                      N_output_arrays, output_array_types,
-                                      (void **) output_arrays_local);
+  /* add the per_point_status[] array the user-supplied table
+     (create a new empty table if the user didn't supply one) */
+  table_handle = local_interp_handle;
+  if (table_handle < 0)
+  {
+    table_handle = Util_TableCreate (UTIL_TABLE_FLAGS_DEFAULT);
+  }
+  Util_TableSetPointer (table_handle, per_point_status, "per_point_status");
+
+  /* now call the local interpolator for all local points and store
+     the results in the intermediate local output arrays */
+  retval = CCTK_InterpLocalUniform (N_dims, table_handle, param_table_handle,
+                                    origin_local, delta, N_points_local,
+                                    interp_coords_type,
+                                    (const void **) interp_coords_local,
+                                    N_input_arrays, input_array_dims,
+                                    input_array_types, input_arrays,
+                                    N_output_arrays, output_array_types,
+                                    (void **) output_arrays_local);
+
+  /* Now check whether the local interpolator returned per-point
+     status information.
+     If so, and there were per-point errors returned by the local
+     interpolator, then reduce the status values of all points from
+     a requesting processor into a single value (ie. reduce
+     per_point_status[N_points_local] into per_proc_status[nprocs]).
+     The minimum overall status value will be taken. */
+  have_per_point_status = Util_TableGetInt (table_handle, &error_point_status,
+                                            "error_point_status") == 1;
+  if (have_per_point_status && error_point_status)
+  {
+    for (proc = 0; proc < pughGH->nprocs; proc++)
+    {
+      for (point = 0; point < myGH->N_points_to[proc]; point++)
+      {
+        if (per_proc_status[proc] > *per_point_status)
+        {
+          per_proc_status[proc] = *per_point_status;
+        }
+        per_point_status++;
+      }
+    }
+  }
+
+  /* remove the internal options from the user-supplied table */
+  if (table_handle == local_interp_handle)
+  {
+    Util_TableDeleteKey (table_handle, "error_point_status");
+    Util_TableDeleteKey (table_handle, "per_point_status");
+  }
+  else
+  {
+    Util_TableDestroy (table_handle);
+  }
+
+  /* If the local interpolator provided per-point status information
+     then do the global reduction on all per-processor status values now.
+     This will then be taken as the return code to the caller. */
+  if (have_per_point_status)
+  {
+    CACTUS_MPI_ERROR (MPI_Allreduce (per_proc_status, local_interp_status,
+                                     pughGH->nprocs, PUGH_MPI_INT, MPI_MIN,
+                                     pughGH->PUGH_COMM_WORLD));
+    retval = local_interp_status[pughGH->myproc];
+  }
 
-    /* don't need these anymore */
+  /* clean up some intermediate arrays */
+  if (N_points_local > 0)
+  {
     for (i = 0; i < N_dims; i++)
     {
       free (interp_coords_local[i]);
     }
   }
-
-  /* clean up some intermediate arrays */
+  free (per_proc_status);
   free (interp_coords_local);
   free (origin_local);
   free (input_arrays);
@@ -1017,7 +1099,12 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
   free (type_desc);
 #endif /* MPI */
 
-  return (retval);
+  /* store this processor's interpolator status in the parameter table */
+  if (param_table_handle >= 0)
+  {
+    Util_TableSetInt (param_table_handle, retval, "local_interpolator_status");
+  }
+  return (0);
 }