Parallelize AEILocalInterp with OpenMP

This leads to a slight change in behaviour. Currently, AEILocalInterp
traverses the list of points sequentially, and aborts when the first
error is encountered. After parallelisation, there is no fixed order
in which the points are traversed, and if several errors are
encountered, any one of the errors may be returned, not necessarily
the first. I am not aware of any thorn that would or should rely on
such an ordering.

This patch also adds "restrict" and "const" statements that may
improve performance as it gives the compiler more information about
dependencies between pointers.


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

1 files changed, 155 insertions, 138 deletions

diff --git a/src/template.c b/src/template.c
index a22401e..9c012c4 100644
--- a/src/template.c
+++ b/src/template.c
@@ -217,43 +217,43 @@
 
 #ifdef HAVE_OP_I
   static
-    void compute_coeffs_I(FP_XYZ, struct COEFFS_STRUCT *coeffs_I);
+    void compute_coeffs_I(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_I);
 #endif
 #ifdef HAVE_OP_DX
   static
-    void compute_coeffs_dx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dx);
+    void compute_coeffs_dx(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dx);
 #endif
 #ifdef HAVE_OP_DY
   static
-    void compute_coeffs_dy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dy);
+    void compute_coeffs_dy(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dy);
 #endif
 #ifdef HAVE_OP_DZ
   static
-    void compute_coeffs_dz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dz);
+    void compute_coeffs_dz(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dz);
 #endif
 #ifdef HAVE_OP_DXX
   static
-    void compute_coeffs_dxx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxx);
+    void compute_coeffs_dxx(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dxx);
 #endif
 #ifdef HAVE_OP_DXY
   static
-    void compute_coeffs_dxy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxy);
+    void compute_coeffs_dxy(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dxy);
 #endif
 #ifdef HAVE_OP_DXZ
   static
-    void compute_coeffs_dxz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxz);
+    void compute_coeffs_dxz(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dxz);
 #endif
 #ifdef HAVE_OP_DYY
   static
-    void compute_coeffs_dyy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyy);
+    void compute_coeffs_dyy(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dyy);
 #endif
 #ifdef HAVE_OP_DYZ
   static
-    void compute_coeffs_dyz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyz);
+    void compute_coeffs_dyz(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dyz);
 #endif
 #ifdef HAVE_OP_DZZ
   static
-    void compute_coeffs_dzz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dzz);
+    void compute_coeffs_dzz(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dzz);
 #endif
 
 /******************************************************************************/
@@ -335,38 +335,38 @@
 
   @@*/
 int FUNCTION_NAME(/***** coordinate system *****/
-		  const CCTK_REAL coord_origin[],
-		  const CCTK_REAL coord_delta[],
+		  const CCTK_REAL *restrict coord_origin,
+		  const CCTK_REAL *restrict coord_delta,
 		  /***** interpolation points *****/
 		  int N_interp_points,
 		  int interp_coords_type_code,
-		  const void* const interp_coords[],
-		  const CCTK_INT N_boundary_points_to_omit[],
-		  const CCTK_REAL boundary_off_centering_tolerance[],
-		  const CCTK_REAL boundary_extrapolation_tolerance[],
+		  const void *restrict const *restrict interp_coords,
+		  const CCTK_INT *restrict N_boundary_points_to_omit,
+		  const CCTK_REAL *restrict boundary_off_centering_tolerance,
+		  const CCTK_REAL *restrict boundary_extrapolation_tolerance,
 		  /***** input arrays *****/
 		  int N_input_arrays,
-		  const CCTK_INT input_array_offsets[],
-		  const CCTK_INT input_array_strides[],
-		  const CCTK_INT input_array_min_subscripts[],
-		  const CCTK_INT input_array_max_subscripts[],
-		  const CCTK_INT input_array_type_codes[],
-		  const void* const input_arrays[],
+		  const CCTK_INT *restrict input_array_offsets,
+		  const CCTK_INT *restrict input_array_strides,
+		  const CCTK_INT *restrict input_array_min_subscripts,
+		  const CCTK_INT *restrict input_array_max_subscripts,
+		  const CCTK_INT *restrict input_array_type_codes,
+		  const void *restrict const *restrict input_arrays,
 		  /***** output arrays *****/
 		  int N_output_arrays,
-		  const CCTK_INT output_array_type_codes[],
-		  void* const output_arrays[],
+		  const CCTK_INT *restrict output_array_type_codes,
+		  void* const *restrict output_arrays,
 		  /***** operation info *****/
-		  const CCTK_INT operand_indices[],
-		  const CCTK_INT operation_codes[],
+		  const CCTK_INT *restrict operand_indices,
+		  const CCTK_INT *restrict operation_codes,
 		  /***** debugging *****/
 		  int debug, FILE* log_fp,
 		  /***** other return results *****/
 		  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[],
-		  struct Jacobian_info* Jacobian_info)
+		  CCTK_INT *restrict const *restrict molecule_positions,
+		  struct Jacobian_info *restrict Jacobian_info)
 {
 /*
  * ***** Naming conventions *****
@@ -824,6 +824,7 @@ if (debug > 0)
 	fflush(stdout);
 	}
 
+#pragma omp parallel for
     for (pt = 0 ; pt < N_interp_points ; ++pt)
     {
     /* this struct holds a molecule-sized piece of a single */
@@ -871,7 +872,7 @@ if (debug > 0)
      * from the interp_coords[] arrays
      * FIXME: Maybe it would be better (faster) to do this
      *	  with N_DIMS open-coded calls on a function?
-     *	  But then we'd have to have a sentinal value
+     *	  But then we'd have to have a sentinel value
      *	  return for the unknown-type-code error case.
      *	  Yuk! :( :(
      */
@@ -885,16 +886,11 @@ if (debug > 0)
 
 	switch	(interp_coords_type_code)
 		{
-	case CCTK_VARIABLE_REAL:
-		  {
-		const CCTK_REAL *const interp_coords_ptr_real
-			= (const CCTK_REAL *) interp_coords_ptr;
-		interp_coords_fp[axis] = interp_coords_ptr_real[pt];
-		break;
-		  }
-
 	#ifdef HAVE_CCTK_REAL4
 	case CCTK_VARIABLE_REAL4:
+        #ifdef CCTK_REAL_PRECISION_4
+	case CCTK_VARIABLE_REAL:
+        #endif
 		  {
 		const CCTK_REAL4 *const interp_coords_ptr_real4
 			= (const CCTK_REAL4 *) interp_coords_ptr;
@@ -905,6 +901,9 @@ if (debug > 0)
 
 	#ifdef HAVE_CCTK_REAL8
 	case CCTK_VARIABLE_REAL8:
+        #ifdef CCTK_REAL_PRECISION_8
+	case CCTK_VARIABLE_REAL:
+        #endif
 		  {
 		const CCTK_REAL8 *const interp_coords_ptr_real8
 			= (const CCTK_REAL8 *) interp_coords_ptr;
@@ -915,6 +914,9 @@ if (debug > 0)
 
 	#ifdef HAVE_CCTK_REAL16
 	case CCTK_VARIABLE_REAL16:
+        #ifdef CCTK_REAL_PRECISION_16
+	case CCTK_VARIABLE_REAL:
+        #endif
 		  {
 		/* FIXME: maybe we should warn (once per cactus run) */
 		/*        that we may be doing arithmetic in lower */
@@ -935,7 +937,9 @@ if (debug > 0)
 			   ,
 			   interp_coords_type_code);
 							/*NOTREACHED*/
-		return UTIL_ERROR_BAD_INPUT;	/*** ERROR RETURN ***/
+#pragma omp critical
+                return_status = UTIL_ERROR_BAD_INPUT;
+                goto next_iteration;
 		/* end of switch (interp_coords_type_code) */
 		}
 
@@ -1152,8 +1156,12 @@ if (debug > 0)
 #endif
 		    }
 
-	    if (this_point_status < return_status)
-	       then return_status = this_point_status;
+	    if (this_point_status < 0)
+	       then {
+#pragma omp critical
+		    if (this_point_status < return_status)
+		       then return_status = this_point_status;
+		    }
 
 	    if (debug >= 6)
 	       then {
@@ -1293,13 +1301,13 @@ if (debug > 0)
 	 * execute the ***load*** the first time in the test at the
 	 * top of the  part  loop below
 	 */
-	const void* input_array_ptr__last_load = NULL;
+	const void *restrict input_array_ptr__last_load = NULL;
 	int part__last_load = -1;
 
 	    for (out = 0 ; out < N_output_arrays ; ++out)
 	    {
 	    const int in = operand_indices[out];
-	    const void* const input_array_ptr = input_arrays[in];
+	    const void *restrict const input_array_ptr = input_arrays[in];
 
 	    /*
  	     * ***decode*** the output array datatype
@@ -1323,7 +1331,9 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 	   ,
 	   out, (int) output_array_type_codes[out], N_output_parts);
 								/*NOTREACHED*/
-		    return UTIL_ERROR_BAD_INPUT;	/*** ERROR RETURN ***/
+#pragma omp critical
+                    return_status = UTIL_ERROR_BAD_INPUT;
+                    goto next_iteration;
 		    }
 
 	      {
@@ -1361,8 +1371,9 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 	   in, (int) input_array_type_codes[in], N_input_parts,
 	   out, (int) output_array_type_codes[out], N_output_parts);
 								/*NOTREACHED*/
-				return UTIL_ERROR_BAD_INPUT;
-							/*** ERROR RETURN ***/
+#pragma omp critical
+                                return_status = UTIL_ERROR_BAD_INPUT;
+                                goto next_iteration;
 				}
 
 			/*
@@ -1375,13 +1386,32 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 					       + input_array_offsets[in];
 			switch	(input_array_type_codes[in])
 				{
+#ifdef HAVE_CCTK_REAL4
+case CCTK_VARIABLE_REAL4:
+#ifdef CCTK_REAL_PRECISION_4
+case CCTK_VARIABLE_REAL:
+#endif
+	  {
+	const CCTK_REAL4 *restrict const input_array_ptr_real4
+		= ((const CCTK_REAL4 *) input_array_ptr) + input_posn;
+	LOAD_DATA_REAL4(input_array_ptr_real4,
+			STRIDE_IJK,
+			&data);
+	break;
+	  }
+#endif
+
+#ifdef HAVE_CCTK_REAL8
+case CCTK_VARIABLE_REAL8:
+#ifdef CCTK_REAL_PRECISION_8
 case CCTK_VARIABLE_REAL:
+#endif
 	  {
-	const CCTK_REAL *const input_array_ptr_real
-		= ((const CCTK_REAL *) input_array_ptr) + input_posn;
-	LOAD_DATA_REAL(input_array_ptr_real,
-		       STRIDE_IJK,
-		       &data);
+	const CCTK_REAL8 *restrict const input_array_ptr_real8
+		= ((const CCTK_REAL8 *) input_array_ptr) + input_posn;
+	LOAD_DATA_REAL8(input_array_ptr_real8,
+                        STRIDE_IJK,
+                        &data);
       #if (N_DIMS == 2) && (MOLECULE_SIZE == 4)
 	if (debug >= 10)
 	   then {
@@ -1408,37 +1438,17 @@ case CCTK_VARIABLE_REAL:
       #endif
 	break;
 	  }
-
-#ifdef HAVE_CCTK_REAL4
-case CCTK_VARIABLE_REAL4:
-	  {
-	const CCTK_REAL4 *const input_array_ptr_real4
-		= ((const CCTK_REAL4 *) input_array_ptr) + input_posn;
-	LOAD_DATA_REAL4(input_array_ptr_real4,
-			STRIDE_IJK,
-			&data);
-	break;
-	  }
-#endif
-
-#ifdef HAVE_CCTK_REAL8
-case CCTK_VARIABLE_REAL8:
-	  {
-	const CCTK_REAL8 *const input_array_ptr_real8
-		= ((const CCTK_REAL8 *) input_array_ptr) + input_posn;
-	LOAD_DATA_REAL8(input_array_ptr_real8,
-			STRIDE_IJK,
-			&data);
-	break;
-	  }
 #endif
 
 #ifdef HAVE_CCTK_REAL16
 case CCTK_VARIABLE_REAL16:
+#ifdef CCTK_REAL_PRECISION_16
+case CCTK_VARIABLE_REAL:
+#endif
 	  {
 	/* FIXME: maybe we should warn (once per cactus run) that we may be */
 	/*        doing arithmetic in lower precision than the data type? */
-	const CCTK_REAL16 *const input_array_ptr_real16
+	const CCTK_REAL16 *restrict const input_array_ptr_real16
 		= ((const CCTK_REAL16 *) input_array_ptr) + input_posn;
 	LOAD_DATA_REAL16(input_array_ptr_real16,
 			 STRIDE_IJK,
@@ -1447,21 +1457,13 @@ case CCTK_VARIABLE_REAL16:
 	  }
 #endif
 
-case CCTK_VARIABLE_COMPLEX:
-	  {
-	const CCTK_REAL (*const input_array_ptr_complex)[COMPLEX_N_PARTS]
-		= ((const CCTK_REAL (*)[COMPLEX_N_PARTS]) input_array_ptr)
-		  + input_posn;
-	LOAD_DATA_COMPLEX(input_array_ptr_complex,
-			  STRIDE_IJK, part,
-			  &data);
-	break;
-	  }
-
 #ifdef HAVE_CCTK_COMPLEX8
 case CCTK_VARIABLE_COMPLEX8:
+#ifdef CCTK_COMPLEX_PRECISION_8
+case CCTK_VARIABLE_COMPLEX:
+#endif
 	  {
-	const CCTK_REAL4 (*const input_array_ptr_complex8)[COMPLEX_N_PARTS]
+	const CCTK_REAL4 (*restrict const input_array_ptr_complex8)[COMPLEX_N_PARTS]
 		= ((const CCTK_REAL4 (*)[COMPLEX_N_PARTS]) input_array_ptr)
 		  + input_posn;
 	LOAD_DATA_COMPLEX8(input_array_ptr_complex8,
@@ -1473,8 +1475,11 @@ case CCTK_VARIABLE_COMPLEX8:
 
 #ifdef HAVE_CCTK_COMPLEX16
 case CCTK_VARIABLE_COMPLEX16:
+#ifdef CCTK_COMPLEX_PRECISION_16
+case CCTK_VARIABLE_COMPLEX:
+#endif
 	  {
-	const CCTK_REAL8 (*const input_array_ptr_complex16)[COMPLEX_N_PARTS]
+	const CCTK_REAL8 (*restrict const input_array_ptr_complex16)[COMPLEX_N_PARTS]
 		= ((const CCTK_REAL8 (*)[COMPLEX_N_PARTS]) input_array_ptr)
 		  + input_posn;
 	LOAD_DATA_COMPLEX16(input_array_ptr_complex16,
@@ -1487,10 +1492,13 @@ case CCTK_VARIABLE_COMPLEX16:
 
 #ifdef HAVE_CCTK_COMPLEX32
 case CCTK_VARIABLE_COMPLEX32:
+#ifdef CCTK_COMPLEX_PRECISION_32
+case CCTK_VARIABLE_COMPLEX:
+#endif
 	  {
 	/* FIXME: maybe we should warn (once per cactus run) that we may be */
 	/*        doing arithmetic in lower precision than the data type? */
-	const CCTK_REAL16 (*const input_array_ptr_complex32)[COMPLEX_N_PARTS]
+	const CCTK_REAL16 (*restrict const input_array_ptr_complex32)[COMPLEX_N_PARTS]
 		= ((const CCTK_REAL16 (*)[COMPLEX_N_PARTS]) input_array_ptr)
 		  + input_posn;
 	LOAD_DATA_COMPLEX32(input_array_ptr_complex32,
@@ -1508,7 +1516,9 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 	   ,
 	   (int) input_array_type_codes[in],
 	   in);						/*NOTREACHED*/
-return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
+#pragma omp critical
+return_status = UTIL_ERROR_BAD_INPUT;
+goto next_iteration;
 				/* end of switch (input_array_type_codes[in]) */
 				}
 			  }
@@ -1608,8 +1618,9 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 	   ,
 	   (int) operation_codes[out],
 	   out);
-				return UTIL_ERROR_BAD_INPUT;
-							/*** ERROR RETURN ***/
+#pragma omp critical
+                                return_status = UTIL_ERROR_BAD_INPUT;
+                                goto next_iteration;
 				/* end of switch (operation_codes[out]) */
 				}
 
@@ -1629,29 +1640,13 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 			switch	(output_array_type_codes[out])
 				{
 
-case CCTK_VARIABLE_REAL:
-	  {
-	CCTK_REAL *const output_array_ptr_real
-		= (CCTK_REAL *) output_arrays[out];
-	if (debug >= 10)
-	   then {
-		if ((pt & (pt-1)) == 0) /* pt is 0 or a power of 2 */
-		   then {
-			printf("                 result addr is %p\n",
-			       (void *) &output_array_ptr_real[pt]);
-			printf("                 previous value there was %g\n",
-			       output_array_ptr_real[pt]);
-			fflush(stdout);
-			}
-		}
-	output_array_ptr_real[pt] = (CCTK_REAL) result;
-	break;
-	  }
-
 #ifdef HAVE_CCTK_REAL4
 case CCTK_VARIABLE_REAL4:
+#ifdef CCTK_REAL_PRECISION_4
+case CCTK_VARIABLE_REAL:
+#endif
 	  {
-	CCTK_REAL4 *const output_array_ptr_real4
+	CCTK_REAL4 *restrict const output_array_ptr_real4
 		= (CCTK_REAL4 *) output_arrays[out];
 	output_array_ptr_real4[pt] = (CCTK_REAL4) result;
 	break;
@@ -1660,36 +1655,48 @@ case CCTK_VARIABLE_REAL4:
 
 #ifdef HAVE_CCTK_REAL8
 case CCTK_VARIABLE_REAL8:
+#ifdef CCTK_REAL_PRECISION_8
+case CCTK_VARIABLE_REAL:
+#endif
 	  {
-	CCTK_REAL8 *const output_array_ptr_real8
+	CCTK_REAL8 *restrict const output_array_ptr_real8
 		= (CCTK_REAL8 *) output_arrays[out];
 	output_array_ptr_real8[pt] = (CCTK_REAL8) result;
+	if (debug >= 10)
+	   then {
+		if ((pt & (pt-1)) == 0) /* pt is 0 or a power of 2 */
+		   then {
+			printf("                 result addr is %p\n",
+			       (void *) &output_array_ptr_real8[pt]);
+			printf("                 previous value there was %.17g\n",
+			       (double) output_array_ptr_real8[pt]);
+			fflush(stdout);
+			}
+		}
 	break;
 	  }
 #endif
 
 #ifdef HAVE_CCTK_REAL16
 case CCTK_VARIABLE_REAL16:
+#ifdef CCTK_REAL_PRECISION_16
+case CCTK_VARIABLE_REAL:
+#endif
 	  {
-	CCTK_REAL16 *const output_array_ptr_real16
+	CCTK_REAL16 *restrict const output_array_ptr_real16
 		= (CCTK_REAL16 *) output_arrays[out];
 	output_array_ptr_real16[pt] = (CCTK_REAL16) result;
 	break;
 	  }
 #endif
 
-case CCTK_VARIABLE_COMPLEX:
-	  {
-	CCTK_REAL (*const output_array_ptr_complex)[COMPLEX_N_PARTS]
-		= (CCTK_REAL (*)[COMPLEX_N_PARTS]) output_arrays[out];
-	output_array_ptr_complex[pt][part] = (CCTK_REAL) result;
-	break;
-	  }
-
 #ifdef HAVE_CCTK_COMPLEX8
 case CCTK_VARIABLE_COMPLEX8:
+#ifdef CCTK_COMPLEX_PRECISION_8
+case CCTK_VARIABLE_COMPLEX:
+#endif
 	  {
-	CCTK_REAL4 (*const output_array_ptr_complex8)[COMPLEX_N_PARTS]
+	CCTK_REAL4 (*restrict const output_array_ptr_complex8)[COMPLEX_N_PARTS]
 		= (CCTK_REAL4 (*)[COMPLEX_N_PARTS]) output_arrays[out];
 	output_array_ptr_complex8[pt][part] = (CCTK_REAL4) result;
 	break;
@@ -1698,8 +1705,11 @@ case CCTK_VARIABLE_COMPLEX8:
 
 #ifdef HAVE_CCTK_COMPLEX16
 case CCTK_VARIABLE_COMPLEX16:
+#ifdef CCTK_COMPLEX_PRECISION_16
+case CCTK_VARIABLE_COMPLEX:
+#endif
 	  {
-	CCTK_REAL8 (*const output_array_ptr_complex16)[COMPLEX_N_PARTS]
+	CCTK_REAL8 (*restrict const output_array_ptr_complex16)[COMPLEX_N_PARTS]
 		= (CCTK_REAL8 (*)[COMPLEX_N_PARTS]) output_arrays[out];
 	output_array_ptr_complex16[pt][part] = (CCTK_REAL8) result;
 	break;
@@ -1708,8 +1718,11 @@ case CCTK_VARIABLE_COMPLEX16:
 
 #ifdef HAVE_CCTK_COMPLEX32
 case CCTK_VARIABLE_COMPLEX32:
+#ifdef CCTK_COMPLEX_PRECISION_32
+case CCTK_VARIABLE_COMPLEX:
+#endif
 	  {
-	CCTK_REAL16 (*const output_array_ptr_complex32)[COMPLEX_N_PARTS]
+	CCTK_REAL16 (*restrict const output_array_ptr_complex32)[COMPLEX_N_PARTS]
 		= (CCTK_REAL16 (*)[COMPLEX_N_PARTS]) output_arrays[out];
 	output_array_ptr_complex32[pt][part] = (CCTK_REAL16) result;
 	break;
@@ -1724,7 +1737,9 @@ default:
 		   ,
 		   (int) output_array_type_codes[out],
 		   out);
-	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
+#pragma omp critical
+        return_status = UTIL_ERROR_BAD_INPUT;
+        goto next_iteration;
 				/* end of switch (output type code) */
 				}
 
@@ -1738,7 +1753,7 @@ default:
 		if ( (Jacobian_info != NULL)
 		     && (Jacobian_info->Jacobian_pointer[out] != NULL))
 		   then {
-			CCTK_REAL *const Jacobian_ptr
+			CCTK_REAL *restrict const Jacobian_ptr
 			   = Jacobian_info->Jacobian_pointer[out]
 			     + Jacobian_info->Jacobian_offset[out]
 			     + Jacobian_info->Jacobian_interp_point_stride*pt
@@ -1828,8 +1843,9 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 	   ,
 	   (int) operation_codes[out],
 	   out);
-				return UTIL_ERROR_BAD_INPUT;
-							/*** ERROR RETURN ***/
+#pragma omp critical
+                                return_status = UTIL_ERROR_BAD_INPUT;
+                                goto next_iteration;
 				/* end of switch(operation_codes[out])*/
 				}
 			/* end of Jacobian-query code */
@@ -1845,6 +1861,7 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
       }
       }
 
+  next_iteration:;
     /* end of  for (pt = ...)  loop */
     }
 
@@ -1909,7 +1926,7 @@ return return_status;
 
 #ifdef HAVE_OP_I
   static
-    void compute_coeffs_I(FP_XYZ, struct COEFFS_STRUCT *coeffs_I)
+    void compute_coeffs_I(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_I)
   {
   #include COEFFS_I_COMPUTE_FILE_NAME
   }
@@ -1917,7 +1934,7 @@ return return_status;
 
 #ifdef HAVE_OP_DX
   static
-    void compute_coeffs_dx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dx)
+    void compute_coeffs_dx(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dx)
   {
   #include COEFFS_DX_COMPUTE_FILE_NAME
   }
@@ -1925,7 +1942,7 @@ return return_status;
 
 #ifdef HAVE_OP_DY
   static
-    void compute_coeffs_dy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dy)
+    void compute_coeffs_dy(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dy)
   {
   #include COEFFS_DY_COMPUTE_FILE_NAME
   }
@@ -1933,7 +1950,7 @@ return return_status;
 
 #ifdef HAVE_OP_DZ
   static
-    void compute_coeffs_dz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dz)
+    void compute_coeffs_dz(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dz)
   {
   #include COEFFS_DZ_COMPUTE_FILE_NAME
   }
@@ -1941,7 +1958,7 @@ return return_status;
 
 #ifdef HAVE_OP_DXX
   static
-    void compute_coeffs_dxx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxx)
+    void compute_coeffs_dxx(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dxx)
   {
   #include COEFFS_DXX_COMPUTE_FILE_NAME
   }
@@ -1949,7 +1966,7 @@ return return_status;
 
 #ifdef HAVE_OP_DXY
   static
-    void compute_coeffs_dxy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxy)
+    void compute_coeffs_dxy(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dxy)
   {
   #include COEFFS_DXY_COMPUTE_FILE_NAME
   }
@@ -1957,7 +1974,7 @@ return return_status;
 
 #ifdef HAVE_OP_DXZ
   static
-    void compute_coeffs_dxz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxz)
+    void compute_coeffs_dxz(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dxz)
   {
   #include COEFFS_DXZ_COMPUTE_FILE_NAME
   }
@@ -1965,7 +1982,7 @@ return return_status;
 
 #ifdef HAVE_OP_DYY
   static
-    void compute_coeffs_dyy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyy)
+    void compute_coeffs_dyy(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dyy)
   {
   #include COEFFS_DYY_COMPUTE_FILE_NAME
   }
@@ -1973,7 +1990,7 @@ return return_status;
 
 #ifdef HAVE_OP_DYZ
   static
-    void compute_coeffs_dyz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyz)
+    void compute_coeffs_dyz(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dyz)
   {
   #include COEFFS_DYZ_COMPUTE_FILE_NAME
   }
@@ -1981,7 +1998,7 @@ return return_status;
 
 #ifdef HAVE_OP_DZZ
   static
-    void compute_coeffs_dzz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dzz)
+    void compute_coeffs_dzz(FP_XYZ, struct COEFFS_STRUCT *restrict coeffs_dzz)
   {
   #include COEFFS_DZZ_COMPUTE_FILE_NAME
   }