start implementing Jacobian computation:

query the interpatch interpolator to get its Jacobian
(not fully implemented yet)


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@586 f88db872-0e4f-0410-b76b-b9085cfa78c5

3 files changed, 540 insertions, 59 deletions

diff --git a/src/patch/FIXME b/src/patch/FIXME
new file mode 100644
index 0000000..8765828
--- /dev/null
+++ b/src/patch/FIXME
@@ -0,0 +1,4 @@
+patch_interp::interp_molecule_minmax_mipar() should also query/check
+- molecule family
+- MSS_is_fn_of_gridfn_values and other new queries
+
diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc
index 91af04e..5a3a03a 100644
--- a/src/patch/patch_interp.cc
+++ b/src/patch/patch_interp.cc
@@ -6,6 +6,9 @@
 // patch_interp::~patch_interp
 // patch_interp::[min,max]_ipar_for_gridfn_data
 // patch_interp::interpolate
+// patch_interp::verify_Jacobian_sparsity_pattern_ok
+// patch_interp::interp_molecule_minmax_mipar
+// patch_interp::interp_molecule_posn
 //
 
 #include <stdio.h>
@@ -54,8 +57,8 @@ patch_interp::patch_interp(const patch_edge& my_edge_in,
 	  min_gfn_(my_patch().ghosted_min_gfn()),
 	  max_gfn_(my_patch().ghosted_max_gfn()),
 	  min_iperp_(min_iperp_in), max_iperp_(max_iperp_in),
-	  min_ipar_(min_ipar_for_gridfn_data()),
-	  max_ipar_(max_ipar_for_gridfn_data()),
+	  min_ipar_(compute_min_ipar()),
+	  max_ipar_(compute_max_ipar()),
 	  min_parindex_used_(min_parindex_used_in),
 	  max_parindex_used_(max_parindex_used_in),
 	  interp_par_(interp_par_in),
@@ -65,9 +68,10 @@ patch_interp::patch_interp(const patch_edge& my_edge_in,
 	  // cf. comments in "patch_interp.hh" on gridfn array subscripting
 	  gridfn_coord_origin_(my_edge().par_map().fp_of_int(min_ipar())),
 	  gridfn_coord_delta_ (my_edge().par_map().delta_fp()),
-	  gridfn_type_codes_ (min_gfn_, max_gfn_),
-	  gridfn_data_ptrs_  (min_gfn_, max_gfn_),
-	  result_buffer_ptrs_(min_gfn_, max_gfn_) // no comma
+	  gridfn_type_codes_      (min_gfn_, max_gfn_),
+	  gridfn_data_ptrs_       (min_gfn_, max_gfn_),
+	  interp_data_buffer_ptrs_(min_gfn_, max_gfn_),
+	  Jacobian_ptrs_          (min_gfn_, max_gfn_) // no comma
 {
 int status;
 
@@ -130,7 +134,7 @@ Util_TableDestroy(interp_par_table_handle_);
 // to our constructor, but needs only my_patch_ and my_edge_ members of
 // this object to be valid.
 //
-int patch_interp::min_ipar_for_gridfn_data() const
+int patch_interp::compute_min_ipar() const
 {
 const ghost_zone& min_par_adjacent_ghost_zone
 	= my_patch()
@@ -140,7 +144,7 @@ return min_par_adjacent_ghost_zone.is_symmetry()
        : my_edge().min_ipar_without_corners();
 }
 
-int patch_interp::max_ipar_for_gridfn_data() const
+int patch_interp::compute_max_ipar() const
 {
 const ghost_zone& max_par_adjacent_ghost_zone
 	= my_patch()
@@ -154,14 +158,29 @@ return max_par_adjacent_ghost_zone.is_symmetry()
 
 //
 // This function interpolates the specified range of ghosted gridfns
-// at all the coordinates specified when we were constructed.
-// It stores the interpolated data in  result_buffer(gfn, iperp,parindex) .
-// It calls  error_exit() if the point is out of range or any other
-// interpolator error is detected.
+// at all the coordinates specified when we were constructed.  It stores
+// the interpolated data in  interp_data_buffer(gfn, iperp,parindex) .
+//
+// If Jacobian_buffer != NULL, then ...
+// this function also stores the Jacobian of the interpolated data
+// with respect to this patch's ghosted gridfns,
+//	partial interp_data_buffer(ghosted_gfn, iperp, parindex)
+//	--------------------------------------------------------
+//	 partial ghosted_gridfn(ghosted_gfn, iperp, posn+mipar)
+// in
+//	(*Jacobian_buffer)(iperp, parindex, mipar)
+// where we implicitly assume the Jacobian to be independent of
+// ghosted_gfn, and where
+//	posn = molecule_posn_buffer(iperp, parindex)
+// is the molecule position as returned by  interp_molecule_posn() .
+//
+// This function calls  error_exit() if any interpolation point is out
+// of range or any other error is detected.
 //
 void patch_interp::interpolate(int ghosted_min_gfn_to_interp,
-				 int ghosted_max_gfn_to_interp,
-				 jtutil::array3d<fp>& result_buffer_array)
+			       int ghosted_max_gfn_to_interp,
+			       jtutil::array3d<fp>& interp_data_buffer,
+			       jtutil::array3d<fp>* Jacobian_buffer = NULL)
 	const
 {
 const int N_dims = 1;
@@ -174,6 +193,46 @@ const CCTK_INT *gridfn_type_codes_ptr
 	= & gridfn_type_codes_(ghosted_min_gfn_to_interp);
 
 //
+// we condition all the Jacobian-query stuff on this flag to ensure that:
+// - there is always a valid "first gridfn" for the Jacobian query
+// - we don't have to worry about deleting table entries that we never
+//   set due to there being no iperp values in this interpolation region
+//
+const bool Jacobian_flag
+	= (    (Jacobian_buffer != NULL)
+	    && (N_gridfns != 0)
+	    && (jtutil::how_many_in_range(min_iperp(), max_iperp()) > 0)    );
+
+int status, status1, status2;
+
+//
+// if we're doing a Jacobian query,
+// set Jacobian stride info in parameter table
+//
+if (Jacobian_flag)
+   then {
+	const int Jacobian_interp_point_stride
+		= Jacobian_buffer->subscript_stride_j();
+	status1 = Util_TableSetInt(interp_par_table_handle_,
+				   Jacobian_interp_point_stride,
+				   "Jacobian_interp_point_stride");
+	const CCTK_INT Jacobian_m_strides[N_dims]
+		= { Jacobian_buffer->subscript_stride_k() };
+	status2 = Util_TableSetIntArray(interp_par_table_handle_,
+					N_dims, Jacobian_m_strides,
+					"Jacobian_m_strides");
+
+	if ((status1 < 0) || (status2 < 0))
+	   then error_exit(ERROR_EXIT,
+"***** patch_interp::interpolate():\n"
+"        can't set Jacobian stride info in interpolator parmameter table!\n"
+"        table handle=%d   error status1=%d status2=%d\n"
+,
+				   interp_par_table_handle_, status1, status2);
+								/*NOTREACHED*/
+	}
+
+//
 // do the interpolations at each iperp
 //
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
@@ -181,17 +240,21 @@ const CCTK_INT *gridfn_type_codes_ptr
 	const int min_parindex = min_parindex_used_(iperp);
 	const int max_parindex = max_parindex_used_(iperp);
 
+	//
 	// interpolation-point coordinates
+	//
 	const CCTK_INT N_interp_points
 		= jtutil::how_many_in_range(min_parindex, max_parindex);
 	const fp* const interp_coords_ptr = & interp_par_(iperp, min_parindex);
 	const void* const interp_coords[N_dims]
 		= { static_cast<const void *>(interp_coords_ptr) };
 
-	// pointers to gridfn data to interpolate, buffer for results
-		for (int gfn = ghosted_min_gfn_to_interp ;
-		     gfn <= ghosted_max_gfn_to_interp ;
-		     ++gfn)
+	//
+	// pointers to gridfn data to interpolate, and to result buffer
+	//
+		for (int ghosted_gfn = ghosted_min_gfn_to_interp ;
+		     ghosted_gfn <= ghosted_max_gfn_to_interp ;
+		     ++ghosted_gfn)
 		{
 		// set up data pointer to --> (iperp,min_ipar) gridfn
 		const int start_irho
@@ -201,35 +264,77 @@ const CCTK_INT *gridfn_type_codes_ptr
 		gridfn_data_ptrs_(gfn)
 			= static_cast<const void *>(
 				& my_patch()
-				  .ghosted_gridfn(gfn, start_irho,start_isigma)
+				  .ghosted_gridfn(ghosted_gfn,
+						  start_irho, start_isigma)
 						   );
-		result_buffer_ptrs_(gfn)
+		interp_data_buffer_ptrs_(gfn)
 			= static_cast<void *>(
-				& result_buffer_array(gfn, iperp,min_parindex)
+				& interp_data_buffer(ghosted_gfn,
+						     iperp,min_parindex)
 					     );
 		}
 	const void* const* const gridfn_data
 		= & gridfn_data_ptrs_(ghosted_min_gfn_to_interp);
-	void* const* const result_buffer
-		= & result_buffer_ptrs_(ghosted_min_gfn_to_interp);
-
-	// do the actual interpolation
-	const int status
-		= CCTK_InterpLocalUniform(N_dims,
-					  interp_handle_,
-					  interp_par_table_handle_,
-					  &gridfn_coord_origin_,
-					  &gridfn_coord_delta_,
-					  N_interp_points,
-						  interp_coords_type_code,
-						  interp_coords,
-					  N_gridfns,
-						  &N_gridfn_data_points,
-						  gridfn_type_codes_ptr,
-						  gridfn_data,
-					  N_gridfns,
-						  gridfn_type_codes_ptr,
-						  result_buffer);
+	void* const* const interp_buffer
+		= & interp_data_buffer_ptrs_(ghosted_min_gfn_to_interp);
+
+	//
+	// if we're doing a Jacobian query,
+	// set up the query itself in the parameter table
+	//
+	if (Jacobian_flag)
+	   then {
+		//
+		// since (we assume) the Jacobian to be independent of
+		// ghosted_gfn, we only query ghosted_min_gfn_to_interp;
+		// we pass null pointers for all other ghosted_gfns
+		//
+			for (int ghosted_gfn = ghosted_min_gfn_to_interp ;
+			     ghosted_gfn <= ghosted_max_gfn_to_interp ;
+			     ++ghosted_gfn)
+			{
+			Jacobian_ptrs_(ghosted_gfn) = NULL;
+			}
+		Jacobian_ptrs_(ghosted_min_gfn_to_interp)
+			= static_cast<CCTK_POINTER>(
+				& (*Jacobian_buffer)(iperp, min_parindex, 0)
+						   );
+		status = Util_TableSetPointerArray
+				 (interp_par_table_handle_,
+				  N_gridfns,
+				  & Jacobian_ptrs_(ghosted_min_gfn_to_interp),
+				  "Jacobian_pointer");
+
+		if (status < 0)
+		   then error_exit(ERROR_EXIT,
+"***** patch_interp::interpolate():\n"
+"        can't set Jacobian pointer in interpolator parmameter table!\n"
+"        table handle=%d   error status=%d\n"
+"        at iperp=%d of [%d,%d]!\n"
+,
+				   interp_par_table_handle_, status,
+				   iperp, min_iperp(), max_iperp());
+		}
+
+	//
+	// interpolate (and optionally query the Jacobian)
+	// for all the parindex values at this iperp
+	//
+	status = CCTK_InterpLocalUniform(N_dims,
+					 interp_handle_,
+					 interp_par_table_handle_,
+					 &gridfn_coord_origin_,
+					 &gridfn_coord_delta_,
+					 N_interp_points,
+					    interp_coords_type_code,
+					    interp_coords,
+					 N_gridfns,
+					    &N_gridfn_data_points,
+					    gridfn_type_codes_ptr,
+					    gridfn_data,
+					 N_gridfns,
+					    gridfn_type_codes_ptr,
+					    interp_buffer);
 	if (status < 0)
 	   then error_exit(ERROR_EXIT,
 "***** patch_interp::interpolate():\n"
@@ -238,4 +343,332 @@ const CCTK_INT *gridfn_type_codes_ptr
 			   status, iperp, min_iperp(), max_iperp());
 								/*NOTREACHED*/
 	}
+
+//
+// if we're doing a Jacobian query,
+// delete query entries from the parameter table
+// (so future interpolator calls don't mistakenly redo this query)
+//
+if (Jacobian_flag)
+   then {
+	status  = Util_TableDeleteKey(interp_par_table_handle_,
+				      "Jacobian_pointer");
+	status1 = Util_TableDeleteKey(interp_par_table_handle_,
+				      "Jacobian_interp_point_stride");
+	status2 = Util_TableDeleteKey(interp_par_table_handle_,
+				      "Jacobian_m_strides");
+	if ((status < 0) || (status1 < 0) || (status2 < 0))
+	   then error_exit(ERROR_EXIT,
+"***** patch_interp::interpolat():\n"
+"        can't delete Jacobian query entries\n"
+"        from interpolator parmameter table!\n"
+"        table handle=%d   error status=%d status1=%d status2=%d\n"
+,
+			   interp_par_table_handle_, status, status1, status2);
+								/*NOTREACHED*/
+	}
+}
+
+//******************************************************************************
+
+//
+// This function verifies that our interpolator has a Jacobian sparsity
+// pattern which we grok: at present this means that the molecules are
+// fixed-sized hypercubes, with size/shape independent of the interpolation
+// coordinates and the floating-point values in the input arrays.
+// (N.b. we don't take derivatives in the interpatch interpolation,
+// so it doesn't matter of the molecule size/shape depends on any
+// operation_codes[] values.)
+//
+// FIXME: we should also query/check the molecule family!
+//
+// If the verification is successful, this function is a no-op.  If not,
+// this function does an error_exit() and does not return to its caller.
+//
+void patch_interp::verify_Jacobian_sparsity_pattern_ok()
+	const
+{
+int status, status1, status2, status3;
+
+//
+// do the actual interpolator query
+// ... explicit arguments all set to no-op values
+//     so we only do the query, not any actual interpolation
+//
+const int N_dims = 1;
+const int N_interp_points = 0;
+const int interp_coords_type_code = CCTK_VARIABLE_REAL;
+const int N_input_arrays = 0;
+const CCTK_INT* input_array_dims = NULL;
+const CCTK_INT* input_array_type_codes = NULL;
+const void* const input_arrays = NULL;
+const int N_output_arrays = 0;
+const CCTK_INT output_array_type_codes = NULL;
+void* const output_arrays = NULL;
+
+status = CCTK_InterpLocalUniform(N_dims,
+				 interp_handle_, interp_par_table_handle_,
+				 &gridfn_coord_origin_, &gridfn_coord_delta_,
+				 N_interp_points,
+					 interp_coords_type_code,
+					 interp_coords,
+				 N_input_arrays,
+					 input_array_dims,
+					 input_array_type_codes,
+					 input_arrays,
+				 N_output_arrays,
+					 output_array_type_codes,
+					 output_arrays);
+if (status < 0)
+   then error_exit(ERROR_EXIT,
+"***** patch_interp::verify_Jacobian_sparsity_pattern_ok():\n"
+"        error return %d from interpolator query!\n"
+,
+		   status);					/*NOTREACHED*/
+
+//
+// get Jacobian-sparsity query results from parameter table
+//
+CCTK_INT MSS_is_fn_of_interp_coords, MSS_is_fn_of_which_output;
+status1 = Util_TableGetInt(interp_par_table_handle_,
+			   &MSS_is_fn_of_interp_coords,
+			   "MSS_is_fn_of_interp_coords");
+status2 = Util_TableGetInt(interp_par_table_handle_,
+			   &MSS_is_fn_of_input_array_values,
+			   "MSS_is_fn_of_input_array_values");
+status3 = Util_TableGetInt(interp_par_table_handle_,
+			   &Jacobian_is_fn_of_input_array_values,
+			   "Jacobian_is_fn_of_input_array_values");
+if ((status1 < 0) || (status2 < 0) || (status3 < 0))
+   then error_exit(ERROR_EXIT,
+"***** patch_interp::verify_Jacobian_sparsity_pattern_ok():\n"
+"        can't get Jacobian sparsity info\n"
+"        from interpolator parmameter table!\n"
+"        table handle=%d\n"
+"        error status1=%d status2=%d status3=%d\n"
+,
+		   interp_par_table_handle_,
+		   status1, status2, status3);			/*NOTREACHED*/
+
+//
+// verify that we grok the Jacobian sparsity pattern
+//
+if ( MSS_is_fn_of_interp_coords || MSS_is_fn_of_input_array_values
+     || Jacobian_is_fn_of_input_array_values )
+   then error_exit(ERROR_EXIT,
+"***** patch_interp::verify_Jacobian_sparsity_pattern_ok():\n"
+"        implementation restriction: we only grok Jacobians with\n"
+"        fixed-sized hypercube-shaped molecules, independent of\n"
+"        the interpolation coordinates and the floating-point values
+"        MSS_is_fn_of_interp_coords=(int)%d\n"
+"        MSS_is_fn_of_input_array_values=(int)%d\n"
+"        Jacobian_is_fn_of_input_array_values=(int)%d\n"
+,
+		   MSS_is_fn_of_interp_coords,
+		   MSS_is_fn_of_input_array_values,
+		   Jacobian_is_fn_of_input_array_values);	/*NOTREACHED*/
+}
+
+//******************************************************************************
+
+//
+// This function queries the interpolator to get the [min,max] m ipar
+// coordinates of the interpolation molecules.
+//
+void patch_interp::interp_molecule_minmax_mipar(int& min_mipar, int& max_mipar)
+	const
+{
+int status, status1, status2;
+
+//
+// set molecule min/max m query entries in parameter table
+//
+status1 = Util_TableSetInt(interp_par_table_handle_,
+			   0, "interp_molecule_min_m");
+status2 = Util_TableSetInt(interp_par_table_handle_,
+			   0, "interp_molecule_max_m");
+if ((status1 < 0) || (status2 < 0) || (status1 < 0) || (status2 < 0))
+   then error_exit(ERROR_EXIT,
+"***** patch_interp::interp_molecule_minmax_mipar():\n"
+"        can't set molecule min/max m queries\n"
+"        in interpolator parmameter table!\n"
+"        table handle=%d   error status1=%d status2=%d\n"
+,
+		   interp_par_table_handle_, status1, status2);	/*NOTREACHED*/
+
+//
+// do the actual interpolator query
+// ... explicit arguments all set to no-op values
+//     so we only do the query, not any actual interpolation
+//
+const int N_dims = 1;
+const int N_interp_points = 0;
+const int interp_coords_type_code = CCTK_VARIABLE_REAL;
+const int N_input_arrays = 0;
+const CCTK_INT* input_array_dims = NULL;
+const CCTK_INT* input_array_type_codes = NULL;
+const void* const input_arrays = NULL;
+const int N_output_arrays = 0;
+const CCTK_INT output_array_type_codes = NULL;
+void* const output_arrays = NULL;
+
+status = CCTK_InterpLocalUniform(N_dims,
+				 interp_handle_, interp_par_table_handle_,
+				 &gridfn_coord_origin_, &gridfn_coord_delta_,
+				 N_interp_points,
+					 interp_coords_type_code,
+					 interp_coords,
+				 N_input_arrays,
+					 input_array_dims,
+					 input_array_type_codes,
+					 input_arrays,
+				 N_output_arrays,
+					 output_array_type_codes,
+					 output_arrays);
+if (status < 0)
+   then error_exit(ERROR_EXIT,
+"***** patch_interp::interp_molecule_minmax_mipar():\n"
+"        error return %d from interpolator query!\n"
+,
+		   status);					/*NOTREACHED*/
+
+//
+// get molecule min/max m query results from parameter table
+// ... must go through temp variables for CCTK_INT to int type conversion
+//
+CCTK_INT interp_molecule_min_m, interp_molecule_max_m;
+status1 = Util_TableGetIntArray(interp_par_table_handle_,
+				N_dims, &interp_molecule_min_m,
+				"interp_molecule_min_m");
+status2 = Util_TableGetIntArray(interp_par_table_handle_,
+				N_dims, &interp_molecule_max_m,
+				"interp_molecule_max_m");
+if ((status1 < 0) || (status2 < 0))
+   then error_exit(ERROR_EXIT,
+"***** patch_interp::interp_molecule_maxmax_mipar():\n"
+"        can't get molecule min/max m query results\n"
+"        from interpolator parmameter table!\n"
+"        table handle=%d   error status1=%d status2=%d\n"
+,
+		   interp_par_table_handle_, status1, status2);	/*NOTREACHED*/
+min_mipar = interp_molecule_min_m;
+max_mipar = interp_molecule_min_m;
+
+//
+// delete Jacobian-sparsity-pattern query entries from the parameter table
+// (so future interpolator calls don't mistakenly redo this query)
+//
+status1 = Util_TableDeleteKey(interp_par_table_handle_,
+			      "interp_molecule_min_m");
+status2 = Util_TableDeleteKey(interp_par_table_handle_,
+			      "interp_molecule_max_m");
+if ((status1 < 0) || (status2 < 0))
+   then error_exit(ERROR_EXIT,
+"***** patch_interp::verify_Jacobian_sparsity_pattern_ok():\n"
+"        can't delete molecule min/max m queries\n"
+"        from interpolator parmameter table!\n"
+"        table handle=%d   error status1=%d status2=%d\n"
+,
+		   interp_par_table_handle_, status1, status2);	/*NOTREACHED*/
+}
+
+//******************************************************************************
+
+//
+// This function queries our interpolator to find out the molecule ipar
+// positions (which we implicitly assume to be independent of ghosted_gfn),
+// and stores these in  molecule_posn_buffer(iperp, parindex) .
+//
+void patch_interp::interp_molecule_posn
+	(jtutil::array2d<CCTK_INT>& molecule_posn_buffer)
+	const
+{
+const int N_dims = 1;
+int status;
+
+	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
+	{
+	const int min_parindex = min_parindex_used_(iperp);
+	const int max_parindex = max_parindex_used_(iperp);
+
+	// set up the molecule-position query in the parameter table
+	CCCTK_POINTER molecule_posn_ptrs[N_dims]
+		= {
+		  static_cast<CCTK_POINTER>(
+			& molecule_posn_buffer(iperp, min_parindex)
+					   )
+		  };
+	status = Util_TableSetPointerArray(interp_par_table_handle_,
+					   N_dims, molecule_posn_ptrs,
+					   "interp_molecule_positions");
+	if (status < 0)
+	   then error_exit(ERROR_EXIT,
+"***** patch_interp::interp_molecule_posn():\n"
+"        can't set molecule position query\n"
+"        in interpolator parmameter table at iperp=%d of [%d,%d]!\n"
+"        error status=%d\n"
+,
+			   iperp, min_iperp(), max_iperp(),
+			   status);				/*NOTREACHED*/
+
+	// interpolation-point coordinates
+	const CCTK_INT N_interp_points
+		= jtutil::how_many_in_range(min_parindex, max_parindex);
+	const int interp_coords_type_code = CCTK_VARIABLE_REAL;
+	const fp* const interp_coords_ptr = & interp_par_(iperp, min_parindex);
+	const void* const interp_coords[N_dims]
+		= { static_cast<const void *>(interp_coords_ptr) };
+
+	// query the interpolator to get the molecule positions
+	// for all parindex at this iperp
+	// ... N_{input,output}_arrays set to zero
+	//     so we only do the query, not any actual interpolation
+	const int N_input_arrays = 0;
+	const CCTK_INT* input_array_dims = NULL;
+	const CCTK_INT* input_array_type_codes = NULL;
+	const void* const input_arrays = NULL;
+	const int N_output_arrays = 0;
+	const CCTK_INT output_array_type_codes = NULL;
+	void* const output_arrays = NULL;
+	status = CCTK_InterpLocalUniform(N_dims,
+					 interp_handle_,
+					 interp_par_table_handle_,
+					 &gridfn_coord_origin_,
+					 &gridfn_coord_delta_,
+					 N_interp_points,
+					    interp_coords_type_code,
+					    interp_coords,
+					 N_input_arrays,
+					    input_array_dims,
+					    input_array_type_codes,
+					    input_arrays,
+					 N_output_arrays,
+					    output_array_type_codes,
+					    output_arrays);
+	if (status < 0)
+	   then error_exit(ERROR_EXIT,
+"***** patch_interp::interp_molecule_posn():\n"
+"        error return %d from interpolator query at iperp=%d of [%d,%d]!\n"
+,
+			   status, iperp, min_iperp(), max_iperp());
+								/*NOTREACHED*/
+	}
+
+//
+// delete molecule-position query entry from the parameter table
+// (so future interpolator calls don't mistakenly redo this query)
+//
+if (jtutil::how_many_in_range(min_iperp(), max_iperp()) > 0)
+   then {
+	status = Util_TableDeleteKey(interp_par_table_handle_,
+				     "interp_molecule_positions");
+	if (status < 0)
+	   then error_exit(ERROR_EXIT,
+"***** patch_interp::interp_molecule_posn():\n"
+"        can't delete molecule position query
+"        from interpolator parmameter table!\n"
+"        table handle=%d   error status=%d\n"
+,
+			   interp_par_table_handle_, status);	/*NOTREACHED*/
+	}
 }
diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index c8cf8a0..9e944be 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -58,7 +58,7 @@
 //
 // For example, we might interpolate at the points
 //            ipar ipar ipar ipar ipar ipar ipar ipar ipar
-//              1    2    3    4    5    6    7    8    9    
+//              1    2    3    4    5    6    7    8    9
 // iperp=10           (2a)   (3b)   (4c)
 // iperp=11          (2d)   (3e)  (4f)   (5g)
 // where the (2a) etc are the interpolation points, with 2 etc being
@@ -112,11 +112,6 @@ public:
 	const patch& my_patch() const { return my_patch_; }
 	const patch_edge& my_edge() const { return my_edge_; }
 
-	// min/max (iperp,ipar) of the gridfn data we'll use
-	int min_iperp() const { return min_iperp_; }
-	int max_iperp() const { return max_iperp_; }
-	int min_ipar() const { return min_ipar_; }
-	int max_ipar() const { return max_ipar_; }
 
 	//
 	// ***** main client interface *****
@@ -124,12 +119,44 @@ public:
 public:
 	// interpolate specified range of ghosted gridfns
 	// at all the coordinates specified when we were constructed,
-	// store interpolated data in  result_buffer
-	// ... result buffer indices are (ghosted_gfn, iperp,parindex)
-	//     n.b. iperp is *OUR* iperp coordinate
+	// store interpolated data in
+	//	interp_data_buffer(ghosted_gfn, iperp, parindex)
+	//
+	// if  Jacobian_buffer != NULL , then ...
+	//    also store Jacobian of interpolated data with respect to
+	//    this patch's ghosted gridfns,
+	//	partial interp_data_buffer(ghosted_gfn, iperp, parindex)
+	//	--------------------------------------------------------
+	//	 partial ghosted_gridfn(ghosted_gfn, iperp, posn+mipar)
+	//    in
+	//	(*Jacobian_buffer)(iperp, parindex, mipar)
+	//    where we implicitly assume the Jacobian to be independent
+	//    of ghosted_gfn, and where
+	//	posn = molecule_posn_buffer(iperp, parindex)
+	//    as returned by  interp_molecule_posn()
 	void interpolate(int ghosted_min_gfn_to_interp,
 			 int ghosted_max_gfn_to_interp,
-			 jtutil::array3d<fp>& result_buffer_array)
+			 jtutil::array3d<fp>& interp_data_buffer,
+			 jtutil::array3d<fp>* Jacobian_buffer = NULL);
+		const;
+
+	// verify (no-op if ok, error_exit() if not) that interpolator
+	// has a Jacobian sparsity pattern which we grok: at present this
+	// means molecules are fixed-sized hypercubes, with size/shape
+	// independent of interpolation coordinates and the floating-point
+	// values in the input arrays
+	void verify_Jacobian_sparsity_pattern_ok() const
+
+	// get [min,max] m ipar coordinates of interpolation molecules
+	void interp_molecule_minmax_mipar(int& min_mipar, int& max_mipar)
+		const;
+
+	// get interpolation molecule ipar positions in
+	//  molecule_posn_buffer(iperp, parindex)
+	// ... array type is CCTK_INT so we can pass by reference
+	//     to interpolator
+	void interp_molecule_posn
+		(jtutil::array2d<CCTK_INT>& molecule_posn_buffer)
 		const;
 
 
@@ -137,16 +164,25 @@ public:
 	// ***** internal functions *****
 	//
 private:
-	// range of ipar in this patch from which we may use gridfn data
-	// in interpolation -- depends on adjacent ghost zones' types as
-	// described in comments to constructor (above), but needs only
-	// my_patch_ and my_edge_ members of this object to be valid
-	int min_ipar_for_gridfn_data() const;
-	int max_ipar_for_gridfn_data() const;
+	// [min,max] iperp for interpolation and gridfn data
+	int min_iperp() const { return min_iperp_; }
+	int max_iperp() const { return max_iperp_; }
+
+	// compute [min,max] ipar in this patch from which we may use
+	// gridfn data in interpolation -- depends on adjacent ghost zones'
+	// types as described in comments to constructor (below), but needs
+	// only my_patch_ and my_edge_ members of this object to be valid
+	int compute_min_ipar() const;
+	int compute_max_ipar() const;
+
+	// min/max (iperp,ipar) of the gridfn data we'll use
+	// (uses precomputed values)
+	int min_ipar() const { return min_ipar_; }
+	int max_ipar() const { return max_ipar_; }
 
 
 	//
-	// ***** constructor, destructor *****
+	// ***** constructor, destructor, et al *****
 	//
 public:
 	//
@@ -198,6 +234,10 @@ private:
 	patch_interp(const patch_interp& rhs);
 	patch_interp& operator=(const patch_interp& rhs);
 
+
+	//
+	// ***** data members *****
+	//
 private:
 	const patch& my_patch_;
 	const patch_edge& my_edge_;
@@ -254,9 +294,13 @@ private:
 	// ... index is (gfn)
 	mutable jtutil::array1d<const void*> gridfn_data_ptrs_;
 
-	// --> start of result buffer for interpolation
+	// --> start of interpolation data buffer for each gridfn
 	//     (reset for each iperp)
 	// ... we do *not* own the pointed-to data!
 	// ... index is (gfn)
-	mutable jtutil::array1d<void*> result_buffer_ptrs_;
+	mutable jtutil::array1d<void*> interp_data_buffer_ptrs_;
+
+	// array of pointers to where Jacobian info should be stored
+	// if we do Jacobian queries
+	mutable jtutil::array1d<void*> Jacobian_ptrs_;
 	};