* add support for returning per-point interpolator status

* drop support for returning status info on which side of the grid
  an outside-the-grid point is outside -- this turned out to be awkward
  to implement in combinatino with per-point status, and I believe it
  was only used by my apparent horizon finder to produce slightly
  more informative error messages in a few cases
* convert some commented-out debugging code to #ifdef on various DEBUG symbols


git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/AEILocalInterp/trunk@9 0f49ee68-0e4f-0410-9b9c-b2c123ded7ef

4 files changed, 267 insertions, 133 deletions

diff --git a/src/InterpLocalUniform.c b/src/InterpLocalUniform.c
index d4fafa8..76882d4 100644
--- a/src/InterpLocalUniform.c
+++ b/src/InterpLocalUniform.c
@@ -10,9 +10,10 @@
 ** InterpLocalUniform - main driver routine
 **
 ** check_boundary_tolerances - check boundary tolerances for consistency
+** get_error_point_info - get per-point error reporting info from par table
 ** get_and_decode_molecule_family - get & decode molecule_family from par table
 ** get_molecule_positions - get molecule_positions from parameter table
-** get_Jacobian_query - get Jacobian-query info from parameter table
+** get_Jacobian_info - get Jacobian-query info from parameter table
 ** set_error_info - set error information in parameter table
 ** set_molecule_structure - set molecule structure info in parameter table
 ** set_molecule_min_max_m - set molecule size info in parameter table
@@ -62,6 +63,9 @@
 #include "Lagrange-maximum-degree/all_prototypes.h"
 #include "Hermite/all_prototypes.h"
 
+#undef AEILOCALINTERP_DEBUG	/* define this for verbose debugging output */
+#undef AEILOCALINTERP_DEBUG2	/* define this for even more verbose debugging output */
+
 /* the rcs ID and its dummy function to use it */
 static const char* rcsid = "$Header$";
 CCTK_FILEVERSION(AEITHorns_AEILocalInterp_src_InterpLocalUniform_c)
@@ -125,6 +129,9 @@ static
 	 const CCTK_REAL boundary_off_centering_tolerance[MAX_N_BOUNDARIES],
 	 const CCTK_REAL boundary_extrapolation_tolerance[MAX_N_BOUNDARIES]);
 static
+  int get_error_point_info(int param_table_handle,
+			   struct error_info *p_error_info);
+static
   int get_and_decode_molecule_family
 	(int param_table_handle,
 	 int buffer_size, char molecule_family_string_buffer[],
@@ -139,7 +146,7 @@ static
 			struct Jacobian_info* p_Jacobian_info);
 static
   int set_error_info(int param_table_handle,
-		     struct error_flags* p_error_flags);
+		     struct error_info* p_error_info);
 static
   int set_molecule_structure
 	(int param_table_handle,
@@ -151,12 +158,11 @@ static
 	 const struct molecule_min_max_m_info* p_molecule_min_max_m_info);
 
 static
-  int get_and_check_INT
-	(int param_table_handle, const char name[],
-	 bool mandatory_flag, int default_value,
-	 bool check_range_flag, int min_value, int max_value,
-				const char max_value_string[],
-	 CCTK_INT* p_value);
+  int get_and_check_INT(int param_table_handle, const char name[],
+			bool mandatory_flag, int default_value,
+			bool check_range_flag, int min_value, int max_value,
+			const char max_value_string[],
+			CCTK_INT* p_value);
 static
   int get_INT_array(int param_table_handle, const char name[],
 		    bool default_flag, int default_value,
@@ -1153,12 +1159,22 @@ check_boundary_tolerances(N_boundaries,
 /**************************************/
 
   {
+struct error_info error_info;
+status = get_error_point_info(param_table_handle,
+			      &error_info);
+if (status != 0)
+   then return status;					/*** ERROR RETURN ***/
+
+/**************************************/
+
+  {
 #define MOLECULE_FAMILY_BUFSIZ	(MAX_MOLECULE_FAMILY_STRLEN+1)
 char molecule_family_string[MOLECULE_FAMILY_BUFSIZ];
 enum molecule_family molecule_family;
-status = get_and_decode_molecule_family(param_table_handle,
-			     MOLECULE_FAMILY_BUFSIZ, molecule_family_string,
-			     &molecule_family);
+status = get_and_decode_molecule_family
+		(param_table_handle,
+		 MOLECULE_FAMILY_BUFSIZ, molecule_family_string,
+		 &molecule_family);
 if (status != 0)
    then return status;					/*** ERROR RETURN ***/
 
@@ -1524,8 +1540,11 @@ if (p_interp_fn == NULL_INTERP_FN_PTR)
 	}
 
 /* call the subfunction to actually do the interpolation */
+#ifdef AEILOCALINTERP_DEBUG
+printf("AEILocalInterp::InterpLocalUniform.c: calling interpolator fn (N_interp_points=%d)\n", N_interp_points);
+fflush(stdout);
+#endif
   {
-struct error_flags error_flags;
 struct molecule_structure_flags molecule_structure_flags;
 const int return_code
 	= (*p_interp_fn)(coord_origin, coord_delta,
@@ -1542,33 +1561,34 @@ const int return_code
 			 N_output_arrays,
 			    output_array_type_codes, output_arrays,
 			    operand_indices, operation_codes,
-			 &error_flags,
+			 &error_info,
 			 &molecule_structure_flags,
 			 p_molecule_min_max_m_info,
 			 p_molecule_positions,
 			 p_Jacobian_info);
+#ifdef AEILOCALINTERP_DEBUG
+printf("AEILocalInterp::InterpLocalUniform.c: back from interpolator fn with return_code=%d\n", return_code);
+fflush(stdout);
+#endif
 
 /******************************************************************************/
 
 /*
- * is an interpolation point out of range?
+ * set any further error status in the parameter table
  */
-if (return_code == CCTK_ERROR_INTERP_POINT_OUTSIDE)
+status = set_error_info(param_table_handle,
+			&error_info);
+if (status != 0)
    then {
-	status = set_error_info(param_table_handle,
-				&error_flags);
-	if (status != 0)
+	free(input_array_offsets);
+	free(operand_indices);
+	free(operation_codes);
+	if (p_Jacobian_info != NULL)
 	   then {
-		free(input_array_offsets);
-		free(operand_indices);
-		free(operation_codes);
-		if (p_Jacobian_info != NULL)
-		   then {
-			free(Jacobian_info.Jacobian_offset);
-			free(Jacobian_info.Jacobian_pointer);
-			}
-		return status;				/*** ERROR RETURN ***/
+		free(Jacobian_info.Jacobian_offset);
+		free(Jacobian_info.Jacobian_pointer);
 		}
+	return status;					/*** ERROR RETURN ***/
 	}
 
 /******************************************************************************/
@@ -1625,6 +1645,10 @@ if (p_Jacobian_info != NULL)
 	free(Jacobian_info.Jacobian_pointer);
 	}
 
+#ifdef AEILOCALINTERP_DEBUG
+printf("AEILocalInterp::InterpLocalUniform.c: returning %d\n", return_code);
+fflush(stdout);
+#endif
 return return_code;					/*** NORMAL RETURN ***/
   }
   }
@@ -1641,6 +1665,7 @@ return return_code;					/*** NORMAL RETURN ***/
   }
   }
   }
+  }
 }
 
 /******************************************************************************/
@@ -1687,6 +1712,56 @@ int ibndry;
 /******************************************************************************/
 
 /*
+ * This function tries to get
+ *	CCTK_POINTER per_point_status
+ * from the parameter table, and sets  found_per_point_status  to
+ * indicate whether or not that key was found.
+ *
+ * If  per_point_status  is found, this function casts it to a CCTK_INT*
+ * and stores it in p_error_info->per_point_status.  If not (if there's
+ * no such key in the table), this function sets
+ * p_error_info->per_point_status to a NULL pointer.
+ *
+ * Results:
+ * This function returns 0 for ok, or the (nonzero) return code for error.
+ */
+static
+  int get_error_point_info(int param_table_handle,
+			   struct error_info *p_error_info)
+{
+CCTK_POINTER per_point_status;
+int status;
+
+status = Util_TableGetPointer(param_table_handle,
+			      &per_point_status,
+			      "per_point_status");
+
+if      (status == 1)
+   then {
+	p_error_info->found_per_point_status = true;
+	p_error_info->per_point_status = (CCTK_INT*) per_point_status;
+	}
+else if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+   then {
+	p_error_info->found_per_point_status = false;
+	p_error_info->per_point_status = NULL;
+	}
+else	{
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): bad \"per_point_status\" table entry!\n"
+"                              error status=%d",
+		   status);
+	return status;					/*** ERROR RETURN ***/
+	}
+
+return 0;						/*** NORMAL RETURN ***/
+}
+
+/******************************************************************************/
+
+/*
  * This function gets
  *	const char molecule_family[]
  * from the parameter table, or sets the default value "cube" if this
@@ -1931,46 +2006,33 @@ return 0;						/*** NORMAL RETURN ***/
 /******************************************************************************/
 
 /*
- * This function sets the error-info entries
- *	CCTK_INT error_pt
- *	CCTK_INT error_ibndry
- *	CCTK_INT error_axis
- *	CCTK_INT error_direction
- * in the parameter table.
+ * If p_error_info->found_per_point_status is true, this function
+ * sets the paramater-table entry
+ *	CCTK_INT error_point_status
+ * to the negative of p_error_info->error_count.
  *
  * Results:
  * This function returns 0 for ok, or the (nonzero) return code for error.
  */
 static
   int set_error_info(int param_table_handle,
-		     struct error_flags* p_error_flags)
+		     struct error_info* p_error_info)
 {
-const int status1 = Util_TableSetInt(param_table_handle,
-				     p_error_flags->error_pt,
-				     "error_pt");
-const int status2 = Util_TableSetInt(param_table_handle,
-				     p_error_flags->error_ibndry,
-				     "error_ibndry");
-const int status3 = Util_TableSetInt(param_table_handle,
-				     p_error_flags->error_axis,
-				     "error_axis");
-const int status4 = Util_TableSetInt(param_table_handle,
-				     p_error_flags->error_direction,
-				     "error_direction");
-if ((status1 < 0) || (status2 < 0) || (status3 < 0) || (status4 < 0))
+if (p_error_info->found_per_point_status)
    then {
-	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-		   __LINE__, __FILE__, CCTK_THORNSTRING,
+	const int status = Util_TableSetInt(param_table_handle,
+					    -p_error_info->error_count,
+					    "error_point_status");
+	if (status < 0)
+	   then CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
 "   CCTK_InterpLocalUniform():\n"
-"        error setting error-info table entry/entries!\n"
-"        error status1=%d status2=%d status3=%d status4=%d"
-		   ,
-		   status1, status2, status3, status4);
-	return   (status1 < 0) ? status1		/*** ERROR RETURN ***/
-	       : (status2 < 0) ? status2		/*** ERROR RETURN ***/
-	       : (status3 < 0) ? status3		/*** ERROR RETURN ***/
-	       :		 status4;		/*** ERROR RETURN ***/
+"        error setting \"error_point_status\" in parameter table!"
+"        error status=%d"
+			   ,
+			   status);				
+		return status;				/*** ERROR RETURN ***/
 	}
 
 return 0;						/*** NORMAL RETURN ***/
@@ -2097,12 +2159,11 @@ return 0;						/*** NORMAL RETURN ***/
  * This function returns 0 for ok, or the (nonzero) return code for error.
  */
 static
-  int get_and_check_INT
-	(int param_table_handle, const char name[],
-	 bool mandatory_flag, int default_value,
-	 bool check_range_flag, int min_value, int max_value,
-				const char max_value_string[],
-	 CCTK_INT* p_value)
+  int get_and_check_INT(int param_table_handle, const char name[],
+			bool mandatory_flag, int default_value,
+			bool check_range_flag, int min_value, int max_value,
+			const char max_value_string[],
+			CCTK_INT* p_value)
 {
 CCTK_INT value;
 
@@ -2281,8 +2342,8 @@ return 0;						/*** NORMAL RETURN ***/
  */
 static
   int get_REAL_array(int param_table_handle, const char name[],
-		      CCTK_REAL default_value,
-		      int N, CCTK_REAL buffer[])
+		     CCTK_REAL default_value,
+		     int N, CCTK_REAL buffer[])
 {
 const int status
 	= Util_TableGetRealArray(param_table_handle,
diff --git a/src/InterpLocalUniform.h b/src/InterpLocalUniform.h
index 28b807a..3515ff0 100644
--- a/src/InterpLocalUniform.h
+++ b/src/InterpLocalUniform.h
@@ -139,12 +139,17 @@ enum	molecule_family
 /* number of real "parts" in a complex number */
 #define COMPLEX_N_PARTS 2
 
-struct	error_flags
+struct	error_info
 	{
-	int error_pt;
-	int error_ibndry;
-	int error_axis;
-	int error_direction;
+	/* did we find  error_point_status  in the parameter table? */
+	bool found_per_point_status;
+
+	/* NULL pointer to skip per-point status, or */
+	/* --> array of size N_interp_points to be set to per-point status */
+	CCTK_INT* per_point_status;
+
+	/* count of the number of points in error */
+	CCTK_INT error_count;
 	};
 
 struct	molecule_structure_flags
diff --git a/src/template.c b/src/template.c
index fe16818..6cce81c 100644
--- a/src/template.c
+++ b/src/template.c
@@ -209,6 +209,9 @@
   @version   $Header$
   @@*/
 
+#undef AEILOCALINTERP_DEBUG	/* define this for verbose debugging output */
+#undef AEILOCALINTERP_DEBUG2	/* define this for even more verbose debugging output */
+
 /******************************************************************************/
 
 /*
@@ -279,11 +282,16 @@
 	If you change the arguments for this function, note that you
 	must also change the prototype in "template.h".
 
-  @var		error_flags
-  @vdesc	If we encounter an out-of-range point and this pointer
-		is non-NULL, we store a description of the out-of-range
-		point in the pointed-to structure.
-  @vtype	struct error_flags* error_flags;
+  @var		error_info
+  @vdesc	If error_info->per_point_status is non-NULL,
+		   then we store per-point interpolator status for
+			each point in  error_info->per_point_status[pt] .
+		We set error_info->error_count to the number of
+		points in error.
+		We set the remaining elements of  error_info
+		to a detailed description of the first point in error
+		(if any).
+  @vtype	struct error_info* error_info;
   @vio		out
 
   @var		molecule_structure_flags
@@ -349,7 +357,7 @@ int FUNCTION_NAME(/***** coordinate system *****/
 		  const CCTK_INT operand_indices[],
 		  const CCTK_INT operation_codes[],
 		  /***** other return results *****/
-		  struct error_flags* error_flags,
+		  struct error_info* error_info,
 		  struct molecule_structure_flags* molecule_structure_flags,
 		  struct molecule_min_max_m_info* molecule_min_max_m_info,
 		  CCTK_INT* const molecule_positions[],
@@ -392,6 +400,9 @@ int FUNCTION_NAME(/***** coordinate system *****/
  *	    }
  *	compute position of interpolation molecules
  *	   with respect to the grid into the XYZ variables
+ *	   and store per-point status if requested
+ *	if (this point is outside the grid)
+ *	   then continue;
  *	if (querying molecule positions)
  *	   then store this molecule position
  *	compute molecule center 1-d position in input arrays
@@ -744,6 +755,12 @@ int out;
  * interpolate at each point
  */
 int pt;
+int return_status = 0;
+error_info->error_count = 0;
+#ifdef AEILOCALINTERP_DEBUG
+printf("AEILocalInterp:: in interpolator fn: N_DIMS=%d N_interp_points=%d\n", (int) N_DIMS, N_interp_points);
+fflush(stdout);
+#endif
     for (pt = 0 ; pt < N_interp_points ; ++pt)
     {
     /* this struct holds a molecule-sized piece of a single */
@@ -859,23 +876,28 @@ int pt;
 		/* end of switch (interp_coords_type_code) */
 		}
 
-	/* end of for (axis = ...) loop */
+	/* end of for (axis = ...) loop to get interp coords */
 	}
       }
 
-#ifdef PICKLE
-  #if (N_DIMS == 1)
-printf("pt=%d interp coord %.16g\n",
-	pt, (double) interp_coords_fp[0]);
-  #endif
-  #if (N_DIMS == 3)
-printf("pt=%d interp coords %.16g %.16g %.16g\n",
-	pt, (double) interp_coords_fp[0],
-	    (double) interp_coords_fp[1],
-	    (double) interp_coords_fp[2]);
+#ifdef AEILOCALINTERP_DEBUG2
+  #if   (N_DIMS == 1)
+  printf("pt=%d interp coord %.16g\n",
+	 pt, (double) interp_coords_fp[0]);
+  #elif (N_DIMS == 2)
+  printf("pt=%d interp coords %.16g %.16g\n",
+	 pt, (double) interp_coords_fp[0],
+	     (double) interp_coords_fp[1],
+  #elif (N_DIMS == 3)
+  printf("pt=%d interp coords %.16g %.16g %.16g\n",
+	 pt, (double) interp_coords_fp[0],
+	     (double) interp_coords_fp[1],
+	     (double) interp_coords_fp[2]);
+  #else
+    #error "N_DIMS may not be > 3!"
   #endif
-#endif
-
+  fflush(stdout);
+#endif	/* AEILOCALINTERP_DEBUG2 */
 
     /*
      * compute position of interpolation molecules with respect to
@@ -891,6 +913,7 @@ printf("pt=%d interp coords %.16g %.16g %.16g\n",
      *	    are "clean" and thus more likely to be register-allocated
      */
       {
+    int this_point_status = 0;
     int center_ijk_temp[MAX_N_DIMS];
     fp  xyz_temp[MAX_N_DIMS];
     int axis;
@@ -898,41 +921,75 @@ printf("pt=%d interp coords %.16g %.16g %.16g\n",
 	{
 	const int ibndry_min = 2*axis;
 	const int ibndry_max = 2*axis + 1;
-	const int status = AEILocalInterp_molecule_posn
-				(coord_origin[axis], coord_delta[axis],
-				 input_array_min_subscripts[axis]
-				   + N_boundary_points_to_omit[ibndry_min],
-				 input_array_max_subscripts[axis]
-				   - N_boundary_points_to_omit[ibndry_max],
-				 MOLECULE_SIZE,
-				 boundary_off_centering_tolerance[ibndry_min],
-				 boundary_off_centering_tolerance[ibndry_max],
-				 boundary_extrapolation_tolerance[ibndry_min],
-				 boundary_extrapolation_tolerance[ibndry_max],
-				 interp_coords_fp[axis],
-				 &center_ijk_temp[axis], &xyz_temp[axis]);
-	if (status < 0)
-	   then {
-		if (status == MOLECULE_POSN_ERROR_GRID_TINY)
-		   then return CCTK_ERROR_INTERP_GRID_TOO_TINY;
-							/*** ERROR RETURN ***/
-		if (error_flags != NULL)
-		   then {
-			error_flags->error_pt = pt;
-			error_flags->error_axis = axis;
-			if (status == MOLECULE_POSN_ERROR_X_LT_MIN)
-			   then {
-				error_flags->error_ibndry = ibndry_min;
-				error_flags->error_direction = -1;
-				}
-			   else {
-				error_flags->error_ibndry = ibndry_max;
-				error_flags->error_direction = +1;
-				}
-			}
-		return CCTK_ERROR_INTERP_POINT_OUTSIDE;	/*** ERROR RETURN ***/
+	const int mp_status = AEILocalInterp_molecule_posn
+				 (coord_origin[axis], coord_delta[axis],
+				  input_array_min_subscripts[axis]
+				    + N_boundary_points_to_omit[ibndry_min],
+				  input_array_max_subscripts[axis]
+				    - N_boundary_points_to_omit[ibndry_max],
+				  MOLECULE_SIZE,
+				  boundary_off_centering_tolerance[ibndry_min],
+				  boundary_off_centering_tolerance[ibndry_max],
+				  boundary_extrapolation_tolerance[ibndry_min],
+				  boundary_extrapolation_tolerance[ibndry_max],
+				  interp_coords_fp[axis],
+				  &center_ijk_temp[axis], &xyz_temp[axis]);
+
+	/*
+	 * unfortunately, the status codes from AEILocalInterp_molecule_posn()
+	 * are different from the ones we need ==> we have to decode them
+	 * to get the status for this axis
+	 */
+	int this_axis_status;
+	switch	(mp_status)
+		{
+	case 0:
+		this_axis_status = 0;
+		break;
+	case MOLECULE_POSN_ERROR_GRID_TINY:
+		this_axis_status = CCTK_ERROR_INTERP_GRID_TOO_SMALL;
+		break;
+	case MOLECULE_POSN_ERROR_X_LT_MIN:
+	case MOLECULE_POSN_ERROR_X_GT_MAX:
+		this_axis_status = CCTK_ERROR_INTERP_POINT_OUTSIDE;
+		break;
+	default:
+		CCTK_VWarn(BUG_MSG_SEVERITY_LEVEL,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): internal error!\n"
+"        unexpected result %d from AEILocalInterp_molecule_posn()\n"
+"        at pt=%d axis=%d!\n"
+			   ,
+			   mp_status, pt, axis);		/*NOTREACHED*/
 		}
+
+	if (this_axis_status < this_point_status)
+	   then this_point_status = this_axis_status;
+
+	/* end of for (axis = ...) loop to locate interp points in the grid */
 	}
+
+    if (error_info->per_point_status != NULL)
+       then error_info->per_point_status[pt] = this_point_status;
+
+    if (this_point_status < 0)
+       then {
+	    ++error_info->error_count;
+
+	    if (this_point_status < return_status)
+	       then return_status = this_point_status;
+
+#ifdef AEILOCALINTERP_DEBUG
+printf("AEILocalInterp:: pt=%d has this_point_status=%d ==> skipping it\n", pt, this_point_status);
+fflush(stdout);
+#endif
+
+	    /* this point is in error (eg outside grid) */
+	    /* ==> don't try to interpolate at this point! */
+	    continue;
+	    }
+
       {
   #if (N_DIMS >= 1)
     const int center_i = center_ijk_temp[X_AXIS];
@@ -947,16 +1004,20 @@ printf("pt=%d interp coords %.16g %.16g %.16g\n",
     const fp  z        = xyz_temp       [Z_AXIS];
   #endif
 
-#ifdef PICKLE
-  #if (N_DIMS == 1)
-printf("interp center_i = %d\n", center_i);
-printf("interp grid-relative x = %.16g\n", (double) x);
-  #endif
-  #if (N_DIMS == 3)
-printf("interp center_ijk = %d %d %d\n", center_i, center_j, center_k);
-printf("interp grid-relative xyz = %.16g %.16g %.16g\n",
-	(double) x, (double) y, (double) z);
+#ifdef AEILOCALINTERP_DEBUG
+  #if   (N_DIMS == 1)
+  printf("interp center_i = %d\n", center_i);
+  printf("interp grid-relative x = %.16g\n", (double) x);
+  #elif (N_DIMS == 2)
+  printf("interp center_ij = %d %d\n", center_i, center_j);
+  printf("interp grid-relative xy = %.16g %.16g\n", (double) x, (double) y);
+  #elif (N_DIMS == 3)
+  printf("interp center_ijk = %d %d %d\n", center_i, center_j, center_k);
+  printf("interp grid-relative xyz = %.16g %.16g %.16g\n",
+	 (double) x, (double) y, (double) z);
+    #error "N_DIMS may not be > 3!"
   #endif
+  fflush(stdout);
 #endif
 
     /*
@@ -1352,6 +1413,13 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 			/*
 			 * store result in output array
 			 */
+			#ifdef AEILOCALINTERP_DEBUG
+			if ((pt & (pt-1)) == 0) /* pt is 0 or a power of 2 */
+			   then {
+				printf("AEILocalInterp:: pt=%d out=%d --> storing result %g\n", pt, out, (double) result);
+				fflush(stdout);
+				}
+			#endif
 			switch	(output_array_type_codes[out])
 				{
 
@@ -1563,7 +1631,7 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
     /* end of  for (pt = ...)  loop */
     }
 
-return 0;						/*** NORMAL RETURN ***/
+return return_status;
   }
   }
 }
diff --git a/src/template.h b/src/template.h
index b5eb5b4..d1f9d92 100644
--- a/src/template.h
+++ b/src/template.h
@@ -33,7 +33,7 @@ int FUNCTION_NAME(/***** coordinate system *****/
 		  const CCTK_INT operand_indices[],
 		  const CCTK_INT operation_codes[],
 		  /***** other return results *****/
-		  struct error_flags* error_flags,
+		  struct error_info* error_info,
 		  struct molecule_structure_flags* molecule_structure_flags,
 		  struct molecule_min_max_m_info* molecule_min_max_m_info,
 		  CCTK_INT* const molecule_positions[],