finish (?) Jacobian computation

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

2 files changed, 110 insertions, 141 deletions

diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc
index 3e08add..0784f74 100644
--- a/src/patch/patch_interp.cc
+++ b/src/patch/patch_interp.cc
@@ -6,6 +6,7 @@
 // patch_interp::~patch_interp
 // patch_interp::compute_[min,max]_ipar
 // patch_interp::interpolate
+//
 // patch_interp::query_interpolator
 // patch_interp::verify_Jacobian_sparsity_pattern_ok
 // patch_interp::molecule_minmax_m_ipar
@@ -37,6 +38,8 @@ using jtutil::error_exit;
 #include "ghost_zone.hh"
 
 //*****************************************************************************
+//*****************************************************************************
+//*****************************************************************************
 
 //
 // This function constructs an  patch_interp  object.
@@ -50,8 +53,8 @@ using jtutil::error_exit;
 //
 patch_interp::patch_interp(const patch_edge& my_edge_in,
 			   int min_iperp_in, int max_iperp_in,
-			   const jtutil::array1d<int>& min_parindex_used_in,
-			   const jtutil::array1d<int>& max_parindex_used_in,
+			   const jtutil::array1d<int>& min_parindex_array_in,
+			   const jtutil::array1d<int>& max_parindex_array_in,
 			   const jtutil::array2d<fp>& interp_par_in,
 			   int interp_handle_in, int interp_par_table_handle_in)
 	: my_patch_(my_edge_in.my_patch()),
@@ -61,8 +64,8 @@ patch_interp::patch_interp(const patch_edge& my_edge_in,
 	  min_iperp_(min_iperp_in), max_iperp_(max_iperp_in),
 	  min_ipar_(compute_min_ipar()),
 	  max_ipar_(compute_max_ipar()),
-	  min_parindex_used_(min_parindex_used_in),
-	  max_parindex_used_(max_parindex_used_in),
+	  min_parindex_array_(min_parindex_array_in),
+	  max_parindex_array_(max_parindex_array_in),
 	  interp_par_(interp_par_in),
 	  interp_handle_(interp_handle_in),
 	  interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)),
@@ -161,7 +164,7 @@ 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 does
 // a separate interpolator call for each iperp.)  This function stores
-// the interpolated data in  interp_data_buffer(gfn, iperp,parindex) .
+// the interpolated data in  data_buffer(gfn, iperp,parindex) .
 //
 void patch_interp::interpolate(int ghosted_min_gfn_to_interp,
 			       int ghosted_max_gfn_to_interp,
@@ -175,7 +178,7 @@ const int N_gridfns = jtutil::how_many_in_range(ghosted_min_gfn_to_interp,
 						ghosted_max_gfn_to_interp);
 const CCTK_INT N_gridfn_data_points
 	= jtutil::how_many_in_range(min_ipar(), max_ipar());
-const int interp_coords_type_code = CCTK_VARIABLE_REAL;	// fp == CCTK_REAL
+const int interp_coords_type_code = CCTK_VARIABLE_REAL;
 const CCTK_INT *gridfn_type_codes_ptr
 	= & gridfn_type_codes_(ghosted_min_gfn_to_interp);
 
@@ -184,12 +187,11 @@ const CCTK_INT *gridfn_type_codes_ptr
 //
 	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);
-
 	//
 	// interpolation-point coordinates
 	//
+	const int min_parindex = min_parindex_array_(iperp);
+	const int max_parindex = max_parindex_array_(iperp);
 	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);
@@ -253,67 +255,88 @@ const CCTK_INT *gridfn_type_codes_ptr
 }
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
 
 //
 // This function queries the interpolator with whatever arguments are
-// in the parameter table.  It passes the specified interpolation-point
-// arguments to the interpolator.  It passes no-op values (a single CCTK_REAL
-// input and output array, but NULL input and output array pointers) for
-// the input/output array arguments to the interpolator.  The result is
-// that the interpolator should do only the desired query, not any actual
-// interpolation operation.
+// in the parameter table.  It specifies arguments such that no actual
+// interpolation is done.
+//
+// In particular, the following interpolator arguments are set up properly
+// based on  iperp :
+//	N_dims
+//	operator_handle, param_table_handle,
+//	coord_origin, coord_delta,
+//	N_interp_points, interp_coords_type_code, interp_coords,
+//	input_array_dims	# specifies the correct par size of the
+//				# patch interpolation region at this iperp
+// The following arguments are set to specify a single input array, but
+// with a NULL data pointer so no actual data is used:
+//	N_input_arrays
+//	input_array_type_codes, input_arrays
+// The following arguments are set to specify a single output array, but
+// with a NULL data pointer so no actual data is stored:
+//	N_output_arrays
+//	output_array_type_codes, output_arrays
+//
+// Arguments:
+// iperp = (in) Specifies where in the patch interpolation region the
+//		interpolator query should be done.
 //
 // Results:
-// This function returns the interpolator status code.
+// If the interpolator returns an error status, this function does an
+//  error_exit()  (and does not return to its caller).  Otherwise, this
+// function returns the interpolator status code.
 //
-int patch_interp::query_interpolator(int N_interp_points = 0,
-				     const void *const *interp_coords = NULL)
+int patch_interp::query_interpolator(const char function_name[], int iperp)
 	const
 {
-#ifdef NOT_YET
-const bool interp_points_flag = (iperp != INT_SENTINAL_VALUE);
-
-xxxxx
+int status;
 const int N_dims = 1;
-const int N_interp_points = interp_points_flag
-			    ? jtutil::how_many_in_range(min_parindex,
-							max_parindex)
-			    : 0;
-xxxx
+
+const int min_parindex = min_parindex_array_(iperp);
+const int max_parindex = max_parindex_array_(iperp);
+const int N_interp_points = jtutil::how_many_in_range(min_parindex,
+						      max_parindex);
 const int interp_coords_type_code = CCTK_VARIABLE_REAL;
+const void* const interp_coords[N_dims]
+	= { static_cast<const void *>(& interp_par_(iperp, min_parindex)) };
+
 const int N_input_arrays = 1;
-const CCTK_INT input_array_dims[N_dims] = { 0 };
+const CCTK_INT input_array_dims[N_dims]
+	= { jtutil::how_many_in_range(min_ipar(), max_ipar()) };
 const CCTK_INT input_array_type_codes[N_input_arrays] = { CCTK_VARIABLE_REAL };
 const void* const input_arrays[N_input_arrays] = { NULL };
+
 const int N_output_arrays = 1;
 const CCTK_INT output_array_type_codes[N_output_arrays]
 	= { CCTK_VARIABLE_REAL };
 void* const output_arrays[N_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::query_interpolator():\n"
+"        on behalf of patch_interp::%s()\n"
+"        error return %d from interpolator query at iperp=%d of [%d,%d]!\n"
+,
+		   function_name,
+		   status, iperp, min_iperp(), max_iperp());	/*NOTREACHED*/
 
-	// 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) };
-
-return 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);
-#else
-return 0;
-#endif
+return status;
 }
 
 //******************************************************************************
@@ -335,7 +358,7 @@ return 0;
 void patch_interp::verify_Jacobian_sparsity_pattern_ok()
 	const
 {
-int status, status1, status2, status3;
+int status1, status2, status3;
 
 //
 // query the interpolator
@@ -343,13 +366,7 @@ int status, status1, status2, status3;
 //     Jacobian metainfo is returned by the interpolator on every call,
 //     even without any explicit request
 //
-status = query_interpolator();
-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*/
+query_interpolator("verify_Jacobian_sparsity_pattern_ok");
 
 //
 // get Jacobian-sparsity query results from parameter table
@@ -405,7 +422,8 @@ if ( MSS_is_fn_of_interp_coords || MSS_is_fn_of_input_array_values
 void patch_interp::molecule_minmax_m_ipar(int& min_m_ipar, int& max_m_ipar)
 	const
 {
-int status, status1, status2;
+const int N_dims = 1;
+int status1, status2;
 
 //
 // set molecule min/max m query entries in parameter table
@@ -426,13 +444,7 @@ if ((status1 < 0) || (status2 < 0))
 //
 // query the interpolator
 //
-status = query_interpolator();
-if (status < 0)
-   then error_exit(ERROR_EXIT,
-"***** patch_interp::molecule_minmax_m_ipar():\n"
-"        error return %d from interpolator query!\n"
-,
-		   status);					/*NOTREACHED*/
+query_interpolator("molecule_minmax_m_ipar");
 
 //
 // get molecule min/max m query results from parameter table
@@ -492,11 +504,10 @@ 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);
+	const int min_parindex = min_parindex_array_(iperp);
 
 	// set up the molecule-position query in the parameter table
-	CCCTK_POINTER molecule_posn_ptrs[N_dims]
+	CCTK_POINTER molecule_posn_ptrs[N_dims]
 	    = { static_cast<CCTK_POINTER>(& posn_buffer(iperp, min_parindex)) };
 	status = Util_TableSetPointerArray(interp_par_table_handle_,
 					   N_dims, molecule_posn_ptrs,
@@ -511,24 +522,10 @@ int status;
 			   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 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 (the interpolator will
-	// store them directly in the  posn_buffer(iperp,*)  array
-	status = query_interpolator(N_interp_points, interp_coords);
-	if (status < 0)
-	   then error_exit(ERROR_EXIT,
-"***** patch_interp::molecule_posn():\n"
-"        error return %d from interpolator query at iperp=%d of [%d,%d]!\n"
-,
-			   status, iperp, min_iperp(), max_iperp());
-								/*NOTREACHED*/
+	// store them directly in the  posn_buffer(iperp,*)  array)
+	query_interpolator("molecule_positions", iperp);
 	}
 
 //
@@ -557,9 +554,9 @@ if (jtutil::how_many_in_range(min_iperp(), max_iperp()) > 0)
 // This function queries the interpolator at each iperp to get 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+m_ipar)
+//	partial interpolate() data_buffer(ghosted_gfn, iperp, parindex)
+//	---------------------------------------------------------------
+//	    partial ghosted_gridfn(ghosted_gfn, iperp, posn+m_ipar)
 //
 // This function stores the Jacobian in
 //	Jacobian_buffer(iperp, parindex, m_ipar)
@@ -582,10 +579,6 @@ void patch_interp::Jacobian(jtutil::array3d<fp>& Jacobian_buffer)
 {
 const int N_dims = 1;
 const int N_gridfns = 1;
-const CCTK_INT N_gridfn_data_points
-	= jtutil::how_many_in_range(min_ipar(), max_ipar());
-const CCTK_INT *gridfn_type_codes_ptr
-	= & gridfn_type_codes_(ghosted_min_gfn_to_interp);
 
 int status, status1, status2;
 
@@ -594,12 +587,12 @@ int status, status1, status2;
 // set Jacobian stride info in parameter table
 //
 const int Jacobian_interp_point_stride
-	= Jacobian_buffer->subscript_stride_j();
+	= 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() };	// = { 1 } in practice
+	= { Jacobian_buffer.subscript_stride_k() };	// = { 1 } in practice
 status2 = Util_TableSetIntArray(interp_par_table_handle_,
 				N_dims, Jacobian_m_strides,
 				"Jacobian_m_strides");
@@ -617,17 +610,7 @@ if ((status1 < 0) || (status2 < 0))
 //
 	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);
-
-	//
-	// 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) };
+	const int min_parindex = min_parindex_array_(iperp);
 
 	//
 	// set up the Jacobian query in the parameter table
@@ -651,24 +634,7 @@ if ((status1 < 0) || (status2 < 0))
 	//
 	// query the interpolator to get the Jacobian for this iperp
 	//
-	const CCTK_INT output_array_type_codes[N_gridfns]
-		= { CCTK_VARIABLE_REAL };
-	CCTK_POINTER const output_arrays[N_gridfns] = { NULL };
-#ifdef NOT_YET
-	status = query_interpolator(N_interp_points, interp_coords,
-				    N_gridfns,
-				    output_array_type_codes,
-				    output_arrays);
-#else
-	status = query_interpolator(N_interp_points, interp_coords);
-#endif
-	if (status < 0)
-	   then error_exit(ERROR_EXIT,
-"***** patch_interp::Jacobian():\n"
-"        error return %d from interpolator query at iperp=%d of [%d,%d]!\n"
-,
-			   status, iperp, min_iperp(), max_iperp());
-								/*NOTREACHED*/
+	query_interpolator("Jacobian", iperp);
 	}
 
 //
diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index 4bd402b..92bc127 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -67,8 +67,8 @@
 // region would be described by
 //	[min,max]_iperp_=[10,11]
 //	[min,max]_ipar_=[1,9]
-//	[min,max]_parindex_used_(10)=[2,5]
-//	[min,max]_parindex_used_(11)=[2,6]
+//	[min,max]_parindex_array_(10)=[2,5]
+//	[min,max]_parindex_array_(11)=[2,6]
 //	interp_par_(10,2) = x[a]
 //	interp_par_(10,3) = x[b]
 //	interp_par_(10,4) = x[c]
@@ -120,7 +120,7 @@ public:
 	// interpolate specified range of ghosted gridfns
 	// at all the coordinates specified when we were constructed,
 	// store interpolated data in
-	//	interp_data_buffer(ghosted_gfn, iperp, parindex)
+	//	data_buffer(ghosted_gfn, iperp, parindex)
 	void interpolate(int ghosted_min_gfn_to_interp,
 			 int ghosted_max_gfn_to_interp,
 			 jtutil::array3d<fp>& data_buffer)
@@ -150,9 +150,9 @@ public:
 
 	// get 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+m_ipar)
+	//	partial interpolate() data_buffer(ghosted_gfn, iperp, parindex)
+	//	---------------------------------------------------------------
+	//	    partial ghosted_gridfn(ghosted_gfn, iperp, posn+m_ipar)
 	// store Jacobian in
 	//	Jacobian_buffer(iperp, parindex, m_ipar)
 	// where we implicitly assume the Jacobian to be independent of
@@ -177,17 +177,20 @@ private:
 	int compute_max_ipar() const;
 
 	// min/max (iperp,ipar) of the gridfn data we'll use
-	// (uses precomputed values)
+	// (uses cached values set up by our constructor)
 	int min_ipar() const { return min_ipar_; }
 	int max_ipar() const { return max_ipar_; }
 
 	// query the interpolator via the parameter table
-	// passing the specified arguments to the interpolator,
-	// and no-op values for all the other interpolator arguments
-	// ... return interpolator status code
-	int query_interpolator(int N_interp_points = 0,
-			       const void* const *interp_coords = NULL)
+	// on behalf of the specified function (name used for error msgs only)
+	// ... error_exit() if error return from interpolator,
+	//     otherwise return interpolator status code
+	int query_interpolator(const char function_name[], int iperp)
 		const;
+	// ... use default iperp for queries that don't care
+	int query_interpolator(const char function_name[])
+		const
+		{ return query_interpolator(function_name, min_iperp_); }
 
 
 	//
@@ -200,7 +203,7 @@ public:
 	//		interpolation region is to belong.
 	// [min,max]_iperp_in = The range of iperp for this interpolation
 	//			region
-	// [min,max]_parindex_used_in(iperp)
+	// [min,max]_parindex_array_in(iperp)
 	//	= [min,max] range of parindex actually used at each iperp.
 	//	  We keep references to these arrays, so they should have
 	//	  lifetimes at last as long as that of this object.
@@ -234,8 +237,8 @@ public:
 	//
 	patch_interp(const patch_edge& my_edge_in,
 		     int min_iperp_in, int max_iperp_in,
-		     const jtutil::array1d<int>& min_parindex_used_in,
-		     const jtutil::array1d<int>& max_parindex_used_in,
+		     const jtutil::array1d<int>& min_parindex_array_in,
+		     const jtutil::array1d<int>& max_parindex_array_in,
 		     const jtutil::array2d<fp>& interp_par_in,
 		     int interp_handle_in, int interp_par_table_handle_in);
 	~patch_interp();
@@ -269,8 +272,8 @@ private:
 	// ... these are references to arrays passed in to our constructor
 	//     ==> we do *not* own them!
 	// ... indices are (iperp)
-	const jtutil::array1d<int>& min_parindex_used_;
-	const jtutil::array1d<int>& max_parindex_used_;
+	const jtutil::array1d<int>& min_parindex_array_;
+	const jtutil::array1d<int>& max_parindex_array_;
 
 	// interp_par_(iperp,parindex)
 	//	= Gives the par coordinates at which we will interpolate;