change the way we handle NULL pointers for the input/output/Jacobian

arrays to give more flexibility:
now you can do just a Jacobian query by passing null input/output pointers,


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

1 files changed, 49 insertions, 32 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/template.c b/src/GeneralizedPolynomial-Uniform/template.c
index 6d6c080..fbb136f 100644
--- a/src/GeneralizedPolynomial-Uniform/template.c
+++ b/src/GeneralizedPolynomial-Uniform/template.c
@@ -314,7 +314,7 @@ int FUNCTION_NAME(/***** coordinate system *****/
  *	   then store molecule min/max m
  *	compute "which derivatives are wanted" flags
  *	precompute 1/dx factors
- *		for each interpolation point
+ *		for (int pt = 0 ; pt < N_interp_point ; ++pt)
  *		{
  *		declare all the coefficients
  *		declare all the data-values variables
@@ -324,25 +324,27 @@ int FUNCTION_NAME(/***** coordinate system *****/
  *		if (querying molecule positions)
  *		   then store this molecule position
  *		compute coefficients for all derivatives which are wanted
- *			for each output array
+ *			for (int out = 0 ; out < N_output_arrays ; ++out)
  *			{
- *			if (    (input_arrays[in] == NULL)
- *			     || (output_arrays[out] == NULL)    )
- *			   then continue;	// skip this interpolation?
+ *			const int in = operand_indices[out];
  *			***decode*** the input/output array datatypes
  *			   to determine whether they're real or complex
- *			   (they must both be the same in this regard),
- *			   and define
- *			      int N_parts = data is complex ? 2 : 1;
- *			int part;
- *				for (part = 0 ; part < N_parts ; ++part)
+ *			   (they must both be the same in this regard), then
+ *			const int N_parts = data is complex ? 2 : 1;
+ *				for (int part = 0 ; part < N_parts ; ++part)
  *				{
- *				if (this output array is computed
- *				    using a different input array
- *				    than the previous one || part != 0)
- *				   then ***fetch*** the input array values
+ *				if ( (input_arrays[in] != NULL)
+ *				     && ( (input_arrays[in] != value at
+ *							       last fetch)
+ *					  || (part != value at last fetch) ) )
+ *				   then {
+ *					save input_arrays[in] and part for
+ *					     "previous value" test above
+ *					***fetch*** the molecule-sized piece
+ *						    of  input_arrays[in][part]
+ *						    at this molecule position
  *					            into local "data" variables
- *				  {
+ *					}
  *				if (output_arrays[out] != NULL)
  *				   then {
  *					fp result;
@@ -382,7 +384,6 @@ int FUNCTION_NAME(/***** coordinate system *****/
  *					case ...
  *						}
  *					}
- *				  }
  *				}
  *			}
  *		}
@@ -408,6 +409,7 @@ int FUNCTION_NAME(/***** coordinate system *****/
  * part = 0-origin index selecting real/imaginary part of a complex number
  */
 
+
 /* number of real "parts" in a complex number */
 #define COMPLEX_N_PARTS	2
 
@@ -425,6 +427,7 @@ int FUNCTION_NAME(/***** coordinate system *****/
   #error "MOLECULE_SIZE inconsistent with MOLECULE_{MIN,MAX}_M!"
 #endif
 
+
 /* input array size, strides, and subscripting computation */
 #if (N_DIMS >= 1)
   const int stride_i = input_array_strides[X_AXIS];
@@ -446,6 +449,7 @@ int FUNCTION_NAME(/***** coordinate system *****/
   #error "N_DIMS must be 1, 2, or 3!"
 #endif
 
+
 /* macros used by machine-generated interpolation coefficient expressions */
 /*
  * FIXME: right now this is used as (eg) RATIONAL(1.0,2.0);
@@ -455,9 +459,6 @@ int FUNCTION_NAME(/***** coordinate system *****/
  */
 #define RATIONAL(num,den)	(num/den)
 
-/*
- ***** execution begins here *****
- */
 
 /*
  * Jacobian structure info
@@ -479,6 +480,7 @@ if (molecule_min_max_m_info != NULL)
 		}
 	}
 
+
 /*
  * compute flags specifying which derivatives are wanted
  */
@@ -513,6 +515,7 @@ if (molecule_min_max_m_info != NULL)
 #ifdef HAVE_OP_DZZ
   bool want_dzz = false;
 #endif
+
   {
 int out;
 	for (out = 0 ; out < N_output_arrays ; ++out)
@@ -566,6 +569,7 @@ int out;
 	}
   }
 
+
 /*
  * save origin/delta variables, precompute 1/delta factors
  * ... in theory we could compute only those factors we're going to use,
@@ -603,6 +607,7 @@ int out;
   #endif
 #endif
 
+
 /*
  * interpolate at each point
  */
@@ -820,7 +825,9 @@ int pt;
 	  #error "N_DIMS must be 1, 2, or 3!"
 	#endif
 
-	/* molecule position queries */
+	/*
+	 * molecule position queries
+	 */
 	if (molecule_positions != NULL)
 	   then {
 		#if (N_DIMS >= 1)
@@ -907,15 +914,18 @@ int pt;
 	 */
 	  {
 	int out;
-	const void *input_array_ptr = NULL;
+
+	// next 2 initializers must be invalid values to make sure we
+	// execute the ***fetch*** the first time in the test at the
+	// top of the  part  loop below
+	const void* input_array_ptr__last_fetch = NULL;
+	const int part__last_fetch = -1;
+
 		for (out = 0 ; out < N_output_arrays ; ++out)
 		{
 		const int in = operand_indices[out];
-		if ((input_arrays[in] == NULL) || (output_arrays[out] == NULL))
-		   then continue;		/* skip this interpolation? */
-		  {
-		const int input_center_sub = input_array_offsets[in]
-					     + center_sub;
+		const int input_center_sub
+			= input_array_offsets[in] + center_sub;
 
 		/*
  		 * ***decode*** the input/output array datatypes
@@ -944,17 +954,25 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 
 		  {
 		const int N_parts = N_input_parts;
+		const void* const input_array_ptr = input_arrays[in];
 		int part;
 			for (part = 0 ; part < N_parts ; ++part)
 			{
-			if ( (input_arrays[in] != input_array_ptr)
-			     || (part != 0) )
+			if ( (input_array_ptr != NULL)
+			     &&
+			     ( (input_array_ptr != input_array_ptr__last_fetch)
+			       || (part != part__last_fetch) ) )
 			   then {
+				// remember when we did the following fetch
+				input_array_ptr__last_fetch = input_array_ptr;
+				part__last_fetch = part;
+
 				/*
-				 * ***fetch*** the input array values
-				 *             into local variables
+				 * ***fetch*** the molecule-sized piece
+				 *             of  input_arrays[in][part]
+				 *             at this molecule position
+				 *             into local "data" variables
 				 */
-				input_array_ptr = input_arrays[in];
 				switch	(input_array_type_codes[in])
 					{
 
@@ -1459,7 +1477,6 @@ if ((Jacobian_info != NULL) && (Jacobian_info->Jacobian_pointer[out] != NULL))
 			/* end of  for (part = ...)  loop */
 			}
 		  }
-		  }
 		/* end of  for (out = ...)  loop */
 		}
 	  }