add a test that output_arrays[out] != NULL before doing the interpolation

--> we can now selectively skip interpolation (eg if we just want a Jacobian)
    by passing a null pointer in the appropriate place


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@52 df1f8a13-aa1d-4dd4-9681-27ded5b42416

1 files changed, 205 insertions, 180 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/template.c b/src/GeneralizedPolynomial-Uniform/template.c
index 145532c..ef2d465 100644
--- a/src/GeneralizedPolynomial-Uniform/template.c
+++ b/src/GeneralizedPolynomial-Uniform/template.c
@@ -341,35 +341,41 @@ int FUNCTION_NAME(/***** coordinate system *****/
  *				   then ***fetch*** the input array values
  *					            into local "data" variables
  *				  {
- *				fp result;
- *				switch	(operation_code)
- *					{
- *				case 0:
- *					result = compute the interpolant
- *						 itself as a linear combination
- *						 of data variables & op=0 coeffs
- *					break;
- *				case 1:
- *					result = compute the interpolant
- *						 itself as a linear combination
- *						 of data variables & op=1 coeffs
- *					result *= inverse_dx;
- *					break;
- *				case ...
+ *				if (output_arrays[out] != NULL)
+ *				   then {
+ *					fp result;
+ *					switch	(operation_code)
+ *						{
+ *					case 0:
+ *						result = compute interpolant
+ *							 as a linear combination
+ *							 of data variables
+ *							 and op=0 coeffs
+ *						break;
+ *					case 1:
+ *						result = compute interpolant
+ *							 as a linear combination
+ *							 of data variables
+ *							 and op=1 coeffs
+ *						result *= inverse_dx;
+ *						break;
+ *					case ...
+ *						}
+ *					***store*** result in output array
  *					}
- *				***store*** result in output array
- *				if (querying Jacobian)
+ *				if (querying Jacobian
+ *				    && (Jacobian_pointer[out] != NULL))
  *				   then {
  *					fp factor;
  *					switch	(operation_code)
  *						{
  *					case 0:
- *						store the op=0 Jacobian values
+ *						store op=0 Jacobian values
  *						break;
  *					case 1:
  *						factor = inverse_dx;
- *						store the op=1 Jacobian values
- *						   (multiplying by  factor )
+ *						store factor
+ *						      * op=1 Jacobian values
  *						break;
  *					case ...
  *						}
@@ -915,7 +921,8 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 	   "\n"
 	   "   can't do real input --> complex output\n"
 	   "   or complex input --> real output interpolation!\n"
-	   "   (0-origin) input #%d datatype = %d, output #%d datatype = %d\n"
+	   "   (0-origin) input  #in=%d  datatype = %d\n"
+	   "   (0-origin) output #out=%d datatype = %d\n"
 	   ,
 	   in, (int) input_array_type_codes[in],
 	   out, (int) output_array_type_codes[out]);	/*NOTREACHED*/
@@ -1110,182 +1117,201 @@ case CCTK_VARIABLE_COMPLEX32:
 
 default:
 	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "input datatype %d not supported!",
-		   input_array_type_codes[in]);		/*NOTREACHED*/
+		   "\n"
+		   "   input datatype %d not supported!\n",
+		   "   (0-origin) input #in=%d"
+		   ,
+		   input_array_type_codes[in],
+		   in);						/*NOTREACHED*/
 	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
+					/* end of switch */
+					/*    (input_array_type_codes[in]) */
 					}
+				/* end of ***fetch*** input array values */
 				}
 
-			/*
-			 * compute the interpolant itself
-			 * as a linear combination of the data variables
-			 */
-			  {
-			fp result;
-			switch	(operation_codes[out])
-				{
-			#ifdef HAVE_OP_I
-			  case 0:
-				#include INTERP_I_COMPUTE_FILE_NAME
-				break;
-			#endif
-			#ifdef HAVE_OP_DX
-			  case 1:
-				#include INTERP_DX_COMPUTE_FILE_NAME
-				result *= inverse_delta_x;
-				break;
-			#endif
-			#ifdef HAVE_OP_DY
-			  case 2:
-				#include INTERP_DY_COMPUTE_FILE_NAME
-				result *= inverse_delta_y;
-				break;
-			#endif
-			#ifdef HAVE_OP_DZ
-			  case 3:
-				#include INTERP_DZ_COMPUTE_FILE_NAME
-				result *= inverse_delta_z;
-				break;
-			#endif
-			#ifdef HAVE_OP_DXX
-			  case 11:
-				#include INTERP_DXX_COMPUTE_FILE_NAME
-				result *= inverse_delta_x2;
-				break;
-			#endif
-			#ifdef HAVE_OP_DXY
-			  case 12:
-			  case 21:
-				#include INTERP_DXY_COMPUTE_FILE_NAME
-				result *= inverse_delta_x * inverse_delta_y;
-				break;
-			#endif
-			#ifdef HAVE_OP_DXZ
-			  case 13:
-			  case 31:
-				#include INTERP_DXZ_COMPUTE_FILE_NAME
-				result *= inverse_delta_x * inverse_delta_z;
-				break;
-			#endif
-			#ifdef HAVE_OP_DYY
-			  case 22:
-				#include INTERP_DYY_COMPUTE_FILE_NAME
-				result *= inverse_delta_y2;
-				break;
-			#endif
-			#ifdef HAVE_OP_DYZ
-			  case 23:
-			  case 32:
-				#include INTERP_DYZ_COMPUTE_FILE_NAME
-				result *= inverse_delta_y * inverse_delta_z;
-				break;
-			#endif
-			#ifdef HAVE_OP_DZZ
-			  case 33:
-				#include INTERP_DZZ_COMPUTE_FILE_NAME
-				result *= inverse_delta_z2;
-				break;
-			#endif
-			default:
+if (output_arrays[out] != NULL)
+   then {
+	/*
+	 * compute the interpolant itself
+	 * as a linear combination of the data variables
+	 */
+	fp result;
+	switch	(operation_codes[out])
+		{
+	#ifdef HAVE_OP_I
+	  case 0:
+		#include INTERP_I_COMPUTE_FILE_NAME
+		break;
+	#endif
+	#ifdef HAVE_OP_DX
+	  case 1:
+		#include INTERP_DX_COMPUTE_FILE_NAME
+		result *= inverse_delta_x;
+		break;
+	#endif
+	#ifdef HAVE_OP_DY
+	  case 2:
+		#include INTERP_DY_COMPUTE_FILE_NAME
+		result *= inverse_delta_y;
+		break;
+	#endif
+	#ifdef HAVE_OP_DZ
+	  case 3:
+		#include INTERP_DZ_COMPUTE_FILE_NAME
+		result *= inverse_delta_z;
+		break;
+	#endif
+	#ifdef HAVE_OP_DXX
+	  case 11:
+		#include INTERP_DXX_COMPUTE_FILE_NAME
+		result *= inverse_delta_x2;
+		break;
+	#endif
+	#ifdef HAVE_OP_DXY
+	  case 12:
+	  case 21:
+		#include INTERP_DXY_COMPUTE_FILE_NAME
+		result *= inverse_delta_x * inverse_delta_y;
+		break;
+	#endif
+	#ifdef HAVE_OP_DXZ
+	  case 13:
+	  case 31:
+		#include INTERP_DXZ_COMPUTE_FILE_NAME
+		result *= inverse_delta_x * inverse_delta_z;
+		break;
+	#endif
+	#ifdef HAVE_OP_DYY
+	  case 22:
+		#include INTERP_DYY_COMPUTE_FILE_NAME
+		result *= inverse_delta_y2;
+		break;
+	#endif
+	#ifdef HAVE_OP_DYZ
+	  case 23:
+	  case 32:
+		#include INTERP_DYZ_COMPUTE_FILE_NAME
+		result *= inverse_delta_y * inverse_delta_z;
+		break;
+	#endif
+	#ifdef HAVE_OP_DZZ
+	  case 33:
+		#include INTERP_DZZ_COMPUTE_FILE_NAME
+		result *= inverse_delta_z2;
+		break;
+	#endif
+	  default:
 CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-	   "Generalized interpolation operation_code %d not supported!",
-	   operation_codes[out]);				/*NOTREACHED*/
-				return UTIL_ERROR_BAD_INPUT;
-							/*** ERROR RETURN ***/
-				/* end of switch (operation_codes[out]) */
-				}
+	   "\n"
+	   "   Generalized interpolation operation_code %d not supported!\n"
+	   "   (0-origin) output #out=%d"
+	   ,
+	   operation_codes[out],
+	   out);						/*NOTREACHED*/
+		return UTIL_ERROR_BAD_INPUT;		/*** ERROR RETURN ***/
+		/* end of switch (operation_codes[out]) */
+		}
 
 
-/*
- * ***store*** the result in the output array
- */
-switch	(output_array_type_codes[out])
-	{
+	/*
+	 * ***store*** the result in the output array
+	 */
+	switch	(output_array_type_codes[out])
+		{
 
-case CCTK_VARIABLE_REAL:
-	  {
-	CCTK_REAL *const output_array_ptr_real
-		= (CCTK_REAL *) output_arrays[out];
-	output_array_ptr_real[pt] = result;
-	break;
-	  }
+	case CCTK_VARIABLE_REAL:
+		  {
+		CCTK_REAL *const output_array_ptr_real
+			= (CCTK_REAL *) output_arrays[out];
+		output_array_ptr_real[pt] = result;
+		break;
+		  }
 
-#ifdef HAVE_CCTK_REAL4
-case CCTK_VARIABLE_REAL4:
-	  {
-	CCTK_REAL4 *const output_array_ptr_real4
-		= (CCTK_REAL4 *) output_arrays[out];
-	output_array_ptr_real4[pt] = result;
-	break;
-	  }
-#endif
+	#ifdef HAVE_CCTK_REAL4
+	case CCTK_VARIABLE_REAL4:
+		  {
+		CCTK_REAL4 *const output_array_ptr_real4
+			= (CCTK_REAL4 *) output_arrays[out];
+		output_array_ptr_real4[pt] = result;
+		break;
+		  }
+	#endif
 
-#ifdef HAVE_CCTK_REAL8
-case CCTK_VARIABLE_REAL8:
-	  {
-	CCTK_REAL8 *const output_array_ptr_real8
-		= (CCTK_REAL8 *) output_arrays[out];
-	output_array_ptr_real8[pt] = result;
-	break;
-	  }
-#endif
+	#ifdef HAVE_CCTK_REAL8
+	case CCTK_VARIABLE_REAL8:
+		  {
+		CCTK_REAL8 *const output_array_ptr_real8
+			= (CCTK_REAL8 *) output_arrays[out];
+		output_array_ptr_real8[pt] = result;
+		break;
+		  }
+	#endif
 
-#ifdef HAVE_CCTK_REAL16
-case CCTK_VARIABLE_REAL16:
-	  {
-	CCTK_REAL16 *const output_array_ptr_real16
-		= (CCTK_REAL16 *) output_arrays[out];
-	output_array_ptr_real16[pt] = result;
-	break;
-	  }
-#endif
+	#ifdef HAVE_CCTK_REAL16
+	case CCTK_VARIABLE_REAL16:
+		  {
+		CCTK_REAL16 *const output_array_ptr_real16
+			= (CCTK_REAL16 *) output_arrays[out];
+		output_array_ptr_real16[pt] = 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] = result;
-	break;
-	  }
+	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] = result;
+		break;
+		  }
 
-#ifdef HAVE_CCTK_COMPLEX8
-case CCTK_VARIABLE_COMPLEX8:
-	  {
-	CCTK_REAL4 (*const output_array_ptr_complex8)[COMPLEX_N_PARTS]
-		= (CCTK_REAL4 (*)[COMPLEX_N_PARTS]) output_arrays[out];
-	output_array_ptr_complex8[pt][part] = result;
-	break;
-	  }
-#endif	/* HAVE_CCTK_COMPLEX8 */
+	#ifdef HAVE_CCTK_COMPLEX8
+	case CCTK_VARIABLE_COMPLEX8:
+		  {
+		CCTK_REAL4 (*const output_array_ptr_complex8)[COMPLEX_N_PARTS]
+			= (CCTK_REAL4 (*)[COMPLEX_N_PARTS]) output_arrays[out];
+		output_array_ptr_complex8[pt][part] = result;
+		break;
+		  }
+	#endif	/* HAVE_CCTK_COMPLEX8 */
 
-#ifdef HAVE_CCTK_COMPLEX16
-case CCTK_VARIABLE_COMPLEX16:
-	  {
-	CCTK_REAL8 (*const output_array_ptr_complex16)[COMPLEX_N_PARTS]
-		= (CCTK_REAL8 (*)[COMPLEX_N_PARTS]) output_arrays[out];
-	output_array_ptr_complex16[pt][part] = result;
-	break;
-	  }
-#endif	/* HAVE_CCTK_COMPLEX16 */
+	#ifdef HAVE_CCTK_COMPLEX16
+	case CCTK_VARIABLE_COMPLEX16:
+		  {
+		CCTK_REAL8 (*const output_array_ptr_complex16)[COMPLEX_N_PARTS]
+			= (CCTK_REAL8 (*)[COMPLEX_N_PARTS]) output_arrays[out];
+		output_array_ptr_complex16[pt][part] = result;
+		break;
+		  }
+	#endif	/* HAVE_CCTK_COMPLEX16 */
 
-#ifdef HAVE_CCTK_COMPLEX32
-case CCTK_VARIABLE_COMPLEX32:
-	  {
-	CCTK_REAL16 (*const output_array_ptr_complex32)[COMPLEX_N_PARTS]
-		= (CCTK_REAL16 (*)[COMPLEX_N_PARTS]) output_arrays[out];
-	output_array_ptr_complex32[pt][part] = result;
-	break;
-	  }
-#endif	/* HAVE_CCTK_COMPLEX32 */
+	#ifdef HAVE_CCTK_COMPLEX32
+	case CCTK_VARIABLE_COMPLEX32:
+		  {
+		CCTK_REAL16 (*const output_array_ptr_complex32)[COMPLEX_N_PARTS]
+			= (CCTK_REAL16 (*)[COMPLEX_N_PARTS]) output_arrays[out];
+		output_array_ptr_complex32[pt][part] = result;
+		break;
+		  }
+	#endif	/* HAVE_CCTK_COMPLEX32 */
 
-default:
-	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "output datatype %d not supported",
-		   output_array_type_codes[out]);		/*NOTREACHED*/
-	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
-	/* end of switch (output type code) */
+	default:
+		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+			   "\n"
+			   "   output datatype %d not supported!\n"
+			   "   (0-origin) output #out=%d"
+			   ,
+			   output_array_type_codes[out],
+			   out);				/*NOTREACHED*/
+		return UTIL_ERROR_BAD_INPUT;		/*** ERROR RETURN ***/
+		/* end of switch (output type code) */
+		}
+
+	/* end of if (output_arrays[out] != NULL) */
 	}
 
+
 /*
  * handle querying the Jacobian
  */
@@ -1399,7 +1425,6 @@ if ((Jacobian_info != NULL) && (Jacobian_info->Jacobian_pointer[out] != NULL))
 	  }
 	/* end of Jacobian-query code */
 	}
-			  }
 			/* end of  for (part = ...)  loop */
 			}
 		  }