hould-be-ready-to-compile version to use Cactus local interpolator

with 1-D pointer arrays reset for each iperp
(--> substantial code simplification; tiny efficiency loss)


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

2 files changed, 87 insertions, 84 deletions

diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc
index e21b32f..42a9198 100644
--- a/src/patch/patch_interp.cc
+++ b/src/patch/patch_interp.cc
@@ -70,7 +70,7 @@ patch_frontier::patch_frontier(const patch_edge& my_edge_in,
 		    ? my_edge_in.max_ipar_with_corners()
 		    : my_edge_in.max_ipar_without_corners())
 	  min_parindex_used_(min_parindex_used_in),
-	  N_interp_points_(min_iperp_in, max_iperp_in),
+	  max_parindex_used_(max_parindex_used_in),
 	  interp_par_(interp_par_in),
 	  ghosted_min_gfn_(ghosted_min_gfn_in),
 	  ghosted_max_gfn_(ghosted_max_gfn_in),
@@ -78,9 +78,7 @@ patch_frontier::patch_frontier(const patch_edge& my_edge_in,
 	  interp_par_table_handle_(Util_TableClone(interp_par_table_handle_in)),
 	  gridfn_coord_origin_(my_edge().par_of_ipar(min_ipar())),
 	  gridfn_coord_delta_(my_edge().par_map().delta_fp()),
-	  gridfn_data_ptrs_(min_iperp_in, max_iperp_in,
-			    ghosted_min_gfn_in, ghosted_max_gfn_in),
-	  gridfn_type_codes_(ghosted_min_gfn_in, ghosted_max_gfn_in)
+	  gridfn_type_codes(ghosted_min_gfn_in, ghosted_max_gfn_in)
 {
 int status;
 
@@ -95,26 +93,6 @@ if (status < 0)
 		   interp_par_table_handle_in,
 		   status);					/*NOTREACHED*/
 
-// number of points to interpolate at each iperp
-	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
-	{
-	N_interp_points_(iperp)
-		= jtutil::how_many_in_range(min_parindex_used_in(iperp),
-					    max_parindex_used_in(iperp));
-	}
-
-// set up pointers to start of gridfn data we can use for each iperp
-	for (int gfn = ghosted_min_gfn() ; gfn <= ghosted_max_gfn() ; ++gfn)
-	{
-	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
-	{
-	const int irho   = my_edge().  irho_of_iperp_ipar(iperp, min_ipar());
-	const int isigma = my_edge().isigma_of_iperp_ipar(iperp, min_ipar());
-	gridfn_data_ptrs_(iperp,gfn)
-		= & my_patch().ghosted_gridfn(gfn, irho,isigma);
-	}
-	}
-
 // set up CCTK_VARIABLE_* type codes for gridfns
 	for (int gfn = ghosted_min_gfn() ; gfn <= ghosted_max_gfn() ; ++gfn)
 	{
@@ -137,7 +115,6 @@ if (status < 0)
 ,
 		   interp_par_table_handle_,
 		   status);					/*NOTREACHED*/
-
 }
 
 //*****************************************************************************
@@ -162,25 +139,68 @@ Util_TableDestroy(interp_par_table_handle_);
 void patch_frontier::interpolate(int ghosted_min_gfn, ghosted_max_gfn,
 				 jtutil::array3d<fp>& result_buffer) const
 {
+//
+// We do a separate interpolation for each iperp.
+//
+// The interpolator accesses the gridfn data with the generic
+// subscripting expression
+//	data[offset + i*stride]
+// where i is a 0-origin point index for each interpolation, and where
+// we get to choose the pointer  data  and the  offset  value for each
+// (gridfn,iperp), and the  stride  for each iperp.  Our strategy is to
+// leave  offset  unspecified so it defaults to 0, use a common stride
+// for all iperp (already set up in our constructor), and specify the
+// appropriate data pointers (--> the start of the gridfn data we should
+// use) for each (gfn,iperp).  With this strategy we don't need to modify
+// the parameter table at all after the constructor.
+//
 const int N_dims = 1;
 const int N_gridfns = jtutil::how_many_in_range(ghosted_min_gfn_to_interp,
 						ghosted_max_gfn_to_interp);
+const int interp_coords_type_code = CCTK_VARIABLE_REAL;	// fp == CCTK_REAL
+const CCTK_INT gridfn_data_dims[N_dims]
+	= { jtutil::how_many_in_range(min_ipar(), max_ipar()) };
+const CCTK_INT *gridfn_type_codes_ptr
+	= & gridfn_type_codes_(ghosted_min_gfn_to_interp);
 
+// do the interpolations
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
 	{
 	const int min_parindex = min_parindex_used_(iperp);
-	const CCTK_INT N_interp_points = N_interp_points_(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 *interp_coords_ptr = & interp_par_(iperp, min_parindex);
 	const void* const interp_coords[N_dims]
-		= { static_cast<const void *>(
-			& interp_par_(iperp,min_parindex)
-					     ) };
-	const void* const* input_arrays
-		= static_cast<const void *const *>(
-			& gridfn_data_ptrs(iperp,ghosted_min_gfn_to_interp)
-						  );
-	// FIXME FIXME input_arrays and output_arrays are dubious
-	const CCTK_INT *gridfn_type_codes_ptr
-		= & gridfn_type_codes_(ghosted_min_gfn_to_interp);
+		= { 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)
+		{
+		const int start_irho
+			= my_edge().  irho_of_iperp_ipar(iperp, min_ipar());
+		const int start_isigma
+			= my_edge().isigma_of_iperp_ipar(iperp, min_ipar());
+		gridfn_data_ptrs_(gfn)
+			= & my_patch()
+			    .ghosted_gridfn(gfn, start_irho,start_isigma);
+		result_buffer_ptrs_(gfn)
+			= & result_buffer(gfn, iperp,min_parindex);
+		}
+	const void* const* gridfn_data
+		= static_cast<const void* const*>(
+			& gridfn_data_ptrs_(ghosted_min_gfn_to_interp)
+						 );
+	void* const result_buffer
+		= static_cast<void* const*>(
+			& result_buffer_ptrs_(ghosted_min_gfn_to_interp)
+					   );
+
+	// do the actual interpolation
 	const int status
 		= CCTK_InterpLocalUniform(N_dims,
 					  interp_handle_,
@@ -188,15 +208,15 @@ const int N_gridfns = jtutil::how_many_in_range(ghosted_min_gfn_to_interp,
 					  gridfn_coord_origin_,
 					  gridfn_coord_delta_,
 					  N_interp_points,
-						  CCTK_VARIABLE_REAL,
+						  interp_coords_type_code,
 						  interp_coords,
 					  N_gridfns,
-						  &N_interp_points,
+						  gridfn_data_dims,
 						  gridfn_type_codes_ptr,
-						  input_arrays,
+						  gridfn_data,
 					  N_gridfns,
 						  gridfn_type_codes_ptr,
-						  output_arrays);
+						  result_buffer);
 	if (status < 0)
 	   then error_exit(ERROR_EXIT,
 "***** patch_frontier::interpolate():\n"
diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index 2c92da2..0857860 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -38,13 +38,12 @@
 //
 // In general we interpolate each gridfn at a number of distinct par
 // for each iperp; the integer "parindex" indexes these points.  We
-// attach no particular semantics to parindex, and it need not be 0-origin
-// or have the same range for each iperp.  [In practice, parindex will
-// be the other patch's ipar coordinate.]  We use
-//    parsub = parindex - min_parindex_used_(iperp)
-// as a 0-origin array subscript indexing the points for each interpolation.
-// We assume that the range of parsub for each iperp is roughly similar,
-// so it's ok to use (iperp,parsub) as a 2-D rectangular index space.
+// attach no particular semantics to parindex, and it need not be
+// 0-origin or have the same range for each iperp.  [In practice,
+// parindex will be the other patch's ipar coordinate.]  However,
+// we assume that the range of parindex is roughly similar for each
+// iperp, so it's ok to use (iperp,parindex) as a 2-D rectangular
+// index space.
 //
 // We use the Cactus local interpolator CCTK_InterpLocalUniform()
 // to do the interpolation.  To minimize interpolator overheads, we
@@ -149,16 +148,12 @@ private:
 	// type as described in comments to constructor (above)
 	int min_ipar_, max_ipar_;
 
-	// min parindex at each iperp
-	// ... this is a reference to an array passed in to our constructor
-	//     ==> we do *not* own this!
-	// ... index is (iperp)
+	// [min,max] parindex at each iperp
+	// ... 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_;
-
-	// number of points to interpolate at each iperp
-	// == 1 + maximum parsub
-	// ... index is (iperp)
-	const jtutil::array1d<int> N_interp_points_;
+	const jtutil::array1d<int>& max_parindex_used_;
 
 	// interp_par_(iperp,parindex)
 	//	= Gives the par coordinates at which we will interpolate;
@@ -187,32 +182,20 @@ private:
 	// (par) origin and delta values of the gridfn data
 	fp gridfn_coord_origin_, gridfn_coord_delta_;
 
-	//
-	// We use (the 0-origin) parsub as an interpolator point index
-	// within each interpolation (i.e. at each value of iperp).
-	// The interpolator accesses the gridfn data with the generic
-	// subscripting expression
-	//	data[offset + parsub*stride]
-	// where we get to choose the pointer  data  and the  offset
-	// value for each (gridfn,iperp), and the stride for each iperp.
-	//
-	// Our strategy will be to leave  offset  unspecified so it
-	// defaults to 0, specify a common stride for all (gfn,iperp),
-	// and specify the appropriate data pointer for each (gfn,iperp).
-	// Thus  data  should point to the start of the gridfn data
-	// we should use for each iperp
-	//
-
-	// --> start of gridfn data we should use for each iperp
-	// subscripts are (iperp,gfn) so [gfn] is contiguous at each iperp
-	// ... n.b. we do *not* own this data!
-	jtutil::array2d<fp *> gridfn_data_ptrs_;
-
-	// CCTK_VARIABLE_ type codes for the gridfns
-	// ... all are actually CCTK_REAL, but we need the array for
-	//     the interpolator
-	// ... datatype is CCTK_INT so we can pass this by reference
-	//     to the interpolator
-	// ... index is (ghosted_gfn)
+	// gridfn type codes for interpolator
+	// ... must be CCTK_INT so we can pass by reference to interpolator
+	// ... index is (gfn)
 	jtutil::array1d<CCTK_INT> gridfn_type_codes_;
-	};
+
+	// --> start of gridfn data to use for interpolation
+	//     (reset for each iperp)
+	// ... we do *not* own the pointed-to data!
+	// ... index is (gfn)
+	jtutil::array1d<fp *> gridfn_data_ptrs_;
+
+	// --> start of result buffer for interpolation
+	//     (reset for each iperp)
+	// ... we do *not* own the pointed-to data!
+	// ... index is (gfn)
+	jtutil::array1d<fp *> result_buffer_ptrs_;
+	}