initial checkin of new LocalInterp thorn

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

1 files changed, 1063 insertions, 0 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
new file mode 100644
index 0000000..5ac372c
--- /dev/null
+++ b/doc/documentation.tex
@@ -0,0 +1,1063 @@
+%version   $Header$
+\documentclass{article}
+
+\def\eqref#1{$(\ref{#1})$}
+\def\cf{cf.\hbox{}}
+\def\ie{i.e.\hbox{}}
+\def\eg{e.g.\hbox{}}
+\def\etal{{\it et~al.\/}}
+
+\def\defn#1{{\bf #1}}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\begin{document}
+\title{LocalInterp}
+\author{Jonathan Thornburg, incorporating code from Thomas Radke}
+%
+% We want CVS to expand the Id keyword on the next line, but we don't
+% want TeX to go into math mode to typeset the expansion (because that
+% wouldn't look as nice in the output), so we use the "$ $" construct
+% to get TeX out of math mode again when typesetting the expansion.
+%
+\date{$ $Id$ $}
+\maketitle
+
+\abstract{
+	 This thorn provides processor-local interpolation of
+	 1-D, 2-D, and 3-D data arrays.  It provides several
+	 different interpolators, and supports both the older
+	 \verb|CCTK_InterpLocal()| API, and the newer
+	 \verb|CCTK_InterpLocalUniform()| and
+	 \verb|CCTK_InterpLocalNonUniform()| APIs.
+	 }  
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Introduction}
+
+This thorn provides processor-local interpolation of 1, 2, and 3-D
+arrays.  At present there are 2~interpolators provided (we hope to
+add other interpolators soon):
+\begin{description}
+\item[Uniform Cartesian]
+	This is the local interpolator which used to live
+	in the PUGHInterp thorn.  It was written in C by
+	Thomas Radke in early 2001, drawing on earlier Fortran code
+	by Paul Walker.  It supports the \verb|CCTK_InterpLocal()| API.
+	It provides Lagrange polynomial interpolation for all
+	combinations of~1, 2, and 3~dimensions, and orders~1,
+	2, and~3.  (It would be easy to add additional dimensions
+	and/or orders if desired.)
+\item[Generalized Polynomial]
+	This interpolator was written in C (plus Maple to generate
+	the coefficients) by Jonathan Thornburg in winter 2001-2002.
+	It supports the \verb|CCTK_InterpLocalUniform()| API.
+	It provides Lagrange polynomial interpolation in
+	1~dimension (orders~1-6) and~2 and 3~dimensions (orders~1-4).
+	(Again, it would be easy to add additional dimensions and/or
+	orders.)  It offers a number of options, described below.
+\end{description}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Terminology}
+
+Within Cactus, each interpolator has a character-string name;
+this is mapped to a Cactus \defn{interpolator handle} by
+\verb|CCTK_InterpHandle()|.  For any given interpolator handle,
+there may be a separate interpolator defined for each of the
+interpolation APIs (both the processor-local ones provided by this
+thorn, and the global ones provided by driver-specific interpolator
+thorns such as PUGHInterp).
+
+Terminology for interpolation seems to differ a bit from one author
+to another.  Here we refer to the \defn{point-centered} interpolation
+of a set of \defn{input arrays} (defining data on a uniformly or
+nonuniformly spaced \defn{grid} of \defn{data points}) to a set of
+\defn{interpolation points} (specified by a corresponding set of
+\defn{interpolation coordinates}), with the results being stored
+in a corresponding set of \defn{output arrays}.  Alternatively,
+we may refer to \defn{cell-centered} interpolation, using a grid
+of \defn{data cells} and a set of \defn{interpolation cells}.
+(This latter terminology is common in hydrodynamics interpolators.)
+
+At present all the interpolators do polynomial interpolation, and
+allow the interpolation of multiple input arrays (to the same set of
+interpolation points) in a single interpolator call, using the basic
+algorithm:
+\begin{verbatim}
+for each interpolation point
+{
+choose a finite difference molecule position
+   somewhere near the interpolation point
+        for each input array
+        {
+        compute an interpolating polynomial which approximates
+           the input data at the molecule points
+        output = polynomial(interpolation point)
+        }
+}
+\end{verbatim}
+In the future, we may add other interpolators where the choice of
+molecule is data-dependent (and thus may vary from one input array to
+the next), and/or where the entire input grid is used in each interpolation.
+
+We define the \defn{order} of the interpolation to be the order of
+the fitted polynomial.  That is, in our terminology, linear interpolation
+is order~1, quadratic is order~2, cubic is order~3, etc.  An order~$n$
+interpolator thus has $O(\Delta x^{n+1})$ errors for smooth input data.
+
+However, because the interpolating polynomial generally changes if
+the interpolation point moves from one grid cell to another, unless
+something special is done the interpolating function isn't smooth,
+\ie{} its The 1st~derivative is generically {\em discontinuous\/},
+with $O(\Delta x^n)$ jump discontinuities each time the interpolating
+polynomial changes.  Correspondingly, the interpolation error is
+generically a ``bump function'' which is zero at each grid point
+and rises to a local maximum in each grid cell.  There are interpolation
+algorithms (\eg{} Hermite polynomial and spline interpolation) which
+give better smoothness, but none of the present interpolators implement
+them. :(
+
+As given in the function reference section of the Cactus User's Guide,
+\verb|interp_coords|, \verb|input_arrays|, and \verb|output_arrays| are
+actually all pointers to arrays of \verb|void *| pointers, since we
+support a number of different Cactus data types.  Internally, the
+interpolator casts these \verb|void *| pointers to \verb|CCTK_REAL *|
+or whatever the correct Cactus data types are.  But when describing
+how the interpolator accesses the various arrays, for simplicity we
+usually gloss over this casting, \ie{}~we pretend that \verb|interp_coords|,
+\verb|input_arrays|, and \verb|output_arrays| are pointers to arrays
+of \verb|CCTK_REAL *| pointers.  (This may become clearer once you
+read the next example.)
+
+We use \verb|pt|, \verb|in|, and \verb|out| as generic 0-origin integer
+subscripts into the arrays of interpolation points, input arrays, and
+output arrays respectively.  We use \verb|(i,j,k)| as a generic
+\verb|N_dims|-vector of integer subscripts into the input array
+\verb|input_arrays[in]|.  (Thus \defn{{\tt |(i,j,k)|} space} refers to
+the grid of data points.)  We usually only write array subscripting
+expressions for the 3-D case; the restrictions/generalizations to
+other dimensions should be obvious.
+
+For example, for 3-D interpolation, the (x,y,z) coordinates of a typical
+interpolation point are given by
+\begin{verbatim}
+x = interp_coords[0][pt]
+y = interp_coords[1][pt]
+z = interp_coords[2][pt]
+\end{verbatim}
+(Notice here that as described above, we've implicitly taken
+\verb|interp_coords| to have the C~type
+\verb|const CCTK_REAL *interp_coords[]|, glossing over the casting
+from its actual C~type of \verb|const void *interp_coords[]|.)
+
+We take \verb|axis| to be an integer specifying a dimension,
+\ie{}~0~for~$x$, 1~for~$y$, 2~for~$z$, \dots.
+
+When describing Jacobians and domains of dependence, it's useful to
+introduce the notion of \defn{molecule position}, a nominal reference
+point for the interpolation molecule in \verb|(i,j,k)| coordinate
+space.  (For example, the molecule position might just be the \verb|(i,j,k)|
+coordinates of the molecule's center.)  We also introduce
+\defn{molecule coordinates} \verb|(mi,mj,jk)|, which are just
+\verb|(i,j,k)| coordinates relative to the molecule position.
+We use \verb|m| or as a generic molecule coordinate.  Thus
+(in notation which should be obvious) a generic molecule operation
+can be written
+\begin{equation}
+\verb|output| = \sum_{\tt m} \verb|coeff[posn+m]| \times \verb|input[posn+m]|
+\end{equation}
+Note that there is no requirement that the output be semantically
+located at the position \verb|posn|!  (This may become clearer once
+you look at the example in
+section~\ref{sect-generic-options/Jacobian/fixed-sized-hypercube}.)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Generic Interpolator Options}
+\label{sect-generic-options}
+
+The newer \verb|CCTK_InterpLocalUniform()| and
+\verb|CCTK_InterpLocalNonUniform()| APIs specify a \defn{parameter table}
+(a Cactus key-value table, manipulated by the \verb|Util_Table*()| APIs)
+as a generic mechanism for passing optional input to, and/or returning
+optional results from, the interpolator.  Different interpolators support
+different options; in this section we describe some options which are,
+or will plausibly, be common to multiple interpolators.
+
+Note that except as described in section~\ref{sect-generic-options/caching}
+(``Caching''), all interpolator arguments and parameters apply only to
+the current interpolator call: there is no visible state kept inside
+the interpolator from one call to another.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Interpolation Order}
+
+The uniform Cartesian interpolator encodes the order in the operator
+name, but other interpolators should use a (mandatory) parameter-table
+parameter
+\begin{verbatim}
+const CCTK_INT order;
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Handling of Out-of-Range Interpolation Points}
+
+By default, interpolators should consider it an error if any interpolation
+point lies outside the grid, \ie{} if the ``interpolation'' is actually
+an extrapolation.  (Some interpolators may apply a small ``fuzz'' to
+this test to help avoid floating-point rounding problems.)
+
+Some interpolators may allow this behavior to be changed by the
+optional parameter
+\begin{verbatim}
+const CCTK_REAL out_of_range_tolerance[N_dims];
+\end{verbatim}
+The semantics of this are as follows:  For each \verb|axis|,
+\begin{description}
+\item[\rm If $\hbox{\tt out\_of\_range\_tolerance[axis]} \ge 0.0$,]
+	then an interpolation point is considered to be ``out of range''
+	if and only if the interpolation point is
+	$> \verb|out_of_range_tolerance[axis]|
+	   \times \verb|coord_delta[axis]|$
+	outside the grid in any coordinate.
+\item[\rm If $\hbox{\tt out\_of\_range\_tolerance[axis]} = -1.0$,]
+	then an interpolation point is considered to be ``out of range''
+	if and only if a centered molecule (or more generally,
+	a molecule whose centering is chosen pretending that
+	the grid is of infinite extent), would require data
+	from outside the grid.
+\end{description}
+Other values of \verb|out_of_range_tolerance[axis]| are illegal.
+
+To provide the ``fuzz'' noted above, \verb|out_of_range_tolerance[]|
+should default to having all elements set to a very small positive
+value, say $10^{-14}$.  (This value should really be scaled with the
+floating-point precision used,%%%
+\footnote{%%%
+	 If this scaling is done, a reasonable value for
+	 {\tt out\_of\_range\_tolerance[]} would be around
+	 $100 \epsilon$, where the \defn{machine epsilon}
+	 $\epsilon$ is defined to be the smallest positive
+	 floating-point number such that $1.0 + \epsilon$
+	 compares ``not equal to'' 1.0 in floating-point
+	 arithmetic.  Machine epsilon values can be found
+	 in the standard header file {\tt <float.h>};
+	 for IEEE single and double precision they are
+	 about $1.19{\times}10^{-7}$ and $2.22{\times}10^{-16}$
+	 respectively.
+	 }%%%
+{} but none of the interpolators do this at present.)
+
+If any interpolation points are out of range, interpolators should
+return the error code \verb|CCTK_ERROR_INTERP_POINT_X_RANGE|, and
+report further details about the error by setting the following
+parameter-table entries:
+\begin{verbatim}
+/* which interpolation point is out of range? */
+/* (value is  pt  for out-of-range point) */
+CCTK_INT out_of_range_pt;
+
+/* in which coordinate axis is the point out of range? */
+/* (value is  axis  for out-of-range point) */
+CCTK_INT out_of_range_axis;
+
+/* on which end of the grid is the point out of range? */
+/* (value is -1 for min, +1 for max) */
+CCTK_INT out_of_range_end;
+\end{verbatim}
+
+Note that if multiple points and/or axes are out of range, different
+interpolators may vary in which out-of-range error they report.  That
+is, it is {\em not\/} necessarily the case that the ``first'' such
+error will be the one reported.%%%
+\footnote{%%%
+	 This is to make life simpler for interpolators
+	 which work in parallel over multiple processors.
+	 }%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Molecule Family}
+\label{sect-generic-options/molecule-family}
+
+An interpolator may support different families/shapes of interpolation
+molecules.  Hypercube-shaped molecules are the simplest and most common
+case, but one could also imagine (say) octagon-shaped molecules in 2-D,
+or some generalization of this in higher numbers of dimensions.
+
+The parameter
+\begin{verbatim}
+const char molecule_family[];   /* set with Util_TableSetString() */
+\end{verbatim}
+may be used to query or change the strategy for choosing the molecules.
+
+The semantics are that if this key is present with the value \verb|""|
+(a 0-length ``empty'' string), this queries the default molecule family:
+the interpolator will (re)set the value to a string describing the default
+molecule family (\verb|"cube"| for hypercube-shaped molecules%%%
+\footnote{%%%
+	 Strictly speaking, these should be called
+	 ``parallelipiped-shaped'' molecules, but this
+	 term is so clumsy (and hard to spell!) that
+	 we just call them hypercube-shaped instead.
+	 }%%%
+).
+If this key is present with a value which is a non-empty string, this
+sets the molecule family/shape based on that string.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Non-Contiguous Input Arrays}
+\label{sect-generic-options/non-contiguous-inputs}
+
+Sometimes one of the input ``arrays'' used by the interpolator might not
+be contiguous in memory.  For example, we might want to do 2-D interpolation
+within a plane of a 3-D grid array, but the plane might or might not be
+contiguous in memory.  (Another example would be that the input array
+might be a member of a compact group.)
+
+One way to do this would be to use the hyperslab API to explicitly compute
+the input data on an appropriate hyperslab, then interpolate within that.
+However, some interpolators may support accessing the appropriate hyperslab
+of the input grid directly inside the interpolator.  If this is supported,
+it's probably considerably cheaper than explicitly computing the hyperslab.
+
+If an interpolator supports this, it should use the following (optional)
+parameter-table entries to specify non-contiguous inputs:
+\begin{verbatim}
+const CCTK_INT input_array_offsets[N_input_arrays];
+
+/* the next 3 table entries are shared by all input arrays */
+const CCTK_INT input_array_strides       [N_dims];
+const CCTK_INT input_array_min_subscripts[N_dims];
+const CCTK_INT input_array_max_subscripts[N_dims];
+\end{verbatim}
+
+Then the interpolator would access the input using the generic
+subscripting expression
+\begin{verbatim}
+input_array[in][offset + i*stride_i + j*stride_j + k*stride_k]
+\end{verbatim}
+where
+\begin{verbatim}
+offset = input_array_offsets[in]
+(stride_i,stride_j,stride_k) = input_array_strides[]
+\end{verbatim}
+\verb|(i,j,k)| run from \verb|input_array_min_subscripts[]|
+to \verb|input_array_max_subscripts[]| inclusive
+(n.b.~this is an {\em inclusive\/} range, \ie{}
+$\verb|min| \le \verb|(i,j,k)| \le \verb|max|$).
+
+The defaults are that each input array is contiguous in memory,
+\ie{} \verb|input_array_offsets[]| = 0,
+\verb|stride| determined from \verb|input_array_dims[]|
+              in the usual Fortran manner,
+\verb|input_array_min_subscripts[]| = all 0s, and
+\verb|input_array_max_subscripts[]| = \verb|input_array_dims[]|-1.
+If the stride and max subscript are both specified explicitly, then the
+explicit \verb|input_array_dims[]| argument to
+\verb|CCTK_InterpLocalUniform()| is ignored.
+
+At present all the interpolators take the output arrays to be
+contiguous 1-D arrays.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Derivatives}
+
+Some interpolators can optionally take derivatives as part of the
+interpolation, \ie{} if we view the input data as being samples of
+a smooth function, then instead of estimating values of that function
+at the interpolation points, the interpolator can instead or additionally
+estimate values of various derivatives of that function at the
+interpolation points.  We don't currently implement it, but one could
+also imagine interpolating more general molecule operations such as
+Laplacians, running integrals, etc.
+
+To specify such operations, an interpolator should use the parameter-table
+entries
+\begin{verbatim}
+const CCTK_INT operand_indices[N_output_arrays];
+const CCTK_INT operation_codes[N_output_arrays];
+\end{verbatim}
+The semantics here are that
+\begin{verbatim}
+output_array[out] = op(input_array[operand_indices[out]])
+\end{verbatim}
+where \verb|op| is an operator specified by the \verb|operation_codes[out]|
+value as described below.
+
+Note that \verb|operand_indices[]| doesn't directly name the inputs,
+but rather gives their (0-origin) subscripts in the list of inputs.
+This allows for a more efficient implementation in the (common) case
+where some of the input arrays have many different operations applied
+to them.
+
+Negative \verb|operation_codes[out]| values  are reserved for future
+expansion.  An \verb|operation_codes[out]| value which is $\ge 0$ is
+taken as a decimal integer encoding a coordinate partial derivative:
+each decimal digit means to take the coordinate partial derivative along
+that (1-origin) axis; the order of the digits in a number is ignored.
+For example:
+\begin{flushleft}
+\begin{tabular}{cl}
+\verb|operation_codes[out]|
+	& What it means								\\
+\hline %-----------------------------------------------------------------------
+0	& no derivative, ``just'' interpolate the input array			\\
+1	& interpolate $\partial \big/ \partial x^1$ of the input array		\\
+2	& interpolate $\partial \big/ \partial x^2$ of the input array		\\
+12 or 21
+	& interpolate $\partial^2 \big/ \partial x^1 \partial x^2$
+	  of the input array							\\
+33	& interpolate $\partial^2 \big/ \partial (x^3)^2$ of the input array	%%%\\
+\end{tabular}
+\end{flushleft}
+
+At present we do {\em not\/} have any \verb|#define|s for the
+operation codes in the Cactus header files.  However, you can avoid
+most of the software-engineering problems of having ``magic numbers''
+for the operation codes, by using the macro
+\begin{verbatim}
+#define DERIV(op)   op
+\end{verbatim}
+to mark all such \verb|operation_codes[]| values in your code.
+There's an example of this in section~\ref{sect-example-derivatives}.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Jacobian and Domain of Dependence}
+\label{sect-generic-options/Jacobian}
+
+Sometimes we want to know the Jacobian of the interpolator,
+\begin{equation}
+\frac{\partial \hbox{\tt output\_array[out][pt]}}
+     {\partial \hbox{\tt input\_array[in][(i,j,k)]}}
+							\label{eqn-Jacobian}
+\end{equation}
+and/or ``just'' the set of \verb|(i,j,k)| for which this is nonzero
+(\ie{} the sparsity structure of the Jacobian, equivalently the domain
+of dependence of the output result, or equivalently the interpolation
+molecule size and shape) for a given \verb|out|, \verb|in|, and \verb|pt|.
+
+The complexity of doing this depends (strongly!) on the structure
+of the Jacobian, and in particular on the answers to the following
+questions:
+\begin{itemize}
+\item	What molecule family is being used?%%%
+\footnote{%%%
+	 For something like a spline interpolator the Jacobian
+	 is generally nonzero for all {\tt |(i,j,k)|}, but
+	 exponentially small for most of them; in this case
+	 the Jacobian-query API would probably specify a cutoff
+	 to return an approximate Jacobian with reasonable sparsity.
+	 }%%%
+\item	Does the interpolator always use the same molecule size
+	and shape (possibly depending on the interpolation order,
+	molecule family, or other such parameters), independent
+	of where the interpolation points are in the grid?%%%
+\footnote{%%%
+	 We can always make the interpolation molecules be
+	 the ``same size and shape'' by padding them with
+	 zero entries to bring them all up to the size of
+	 the largest molecule used anywhere in the grid,
+	 but this would be very inefficient if many molecules
+	 were padded in this way.
+	 }%%%
+\item	Does the interpolator use the same interpolation molecule
+	size and shape for each output array?  (The answer to this
+	question may depend on whether and/or what derivatives are
+	being computed.)
+\end{itemize}
+Because the different cases differ so much in their complexity,
+we define several distinct APIs for obtaining the interpolator's
+Jacobian and/or domain of dependence.
+
+%%%%%%%%%%%%%%%%%%%%
+
+\subsubsection{Querying the Interpolator about the Jacobian's structure}
+
+To allow generic code to determine which of the different Jacobian-structure
+cases applies, (and thus which APIs to use), an interpolator which
+supports Jacobian operations should support using the following
+parameter-table entries to query the interpolator about the Jacobian
+structure:
+
+The parameter \verb|molecule_family| may may be used to query what
+molecule family is being used.  This is described in detail in
+section~\ref{sect-generic-options/molecule-family}.
+
+This parameter may be used to query whether the interpolation molecule's
+size and/or shape varies with the interpolation coordinates:
+\begin{verbatim}
+CCTK_INT MSS_is_fn_of_interp_coords;
+\end{verbatim}
+The semantics of this are that if this key is present%%%
+\footnote{%%%
+	 ``MSS'' abbreviates ``{\bf m}olecule {\bf s}ize
+	 and/or {\bf s}hape''.
+	 }%%%
+{} (the value doesn't matter), then the interpolator will (re)set
+the value to 0~if the molecule size and shape do {\em not\/} vary
+with the interpolation coordinates, or 1~if the molecule size and/or
+shape {\em do\/} vary with the interpolation coordinates.
+
+This parameter may be used to query whether the interpolation molecule's
+size and/or shape varies from one output array to another:
+\begin{verbatim}
+CCTK_INT MSS_is_fn_of_which_output;
+\end{verbatim}
+The semantics of this are that if this key is present (the value
+doesn't matter), then the interpolator will (re)set the value to
+0~if the molecule size and shape do {\em not\/} vary from one output
+array to another (this includes the case where there is only 1~output
+array!), or 1~if the molecule size and/or shape {\em do\/} vary from
+one output array to another.  Note that since the answer to this may
+depend on whether and/or what derivatives are being computed, setting
+this key may force the interpolator to scan completely through the
+derivative specifications.  This probably isn't {\em very\/} costly,
+but you may want to avoid unnecessarily paying that cost by deleting
+the query key (\verb|MSS_is_fn_of_which_output|) from the table
+on future interpolator calls.
+
+%%%%%%%%%%%%%%%%%%%%
+
+\subsubsection{Fixed-Size Hypercube-Shaped Molecules}
+\label{sect-generic-options/Jacobian/fixed-sized-hypercube}
+
+The simplest case (and the only one for which we have defined an API
+at present) is when the molecules are hypercube-shaped and of (small)
+fixed size, independent of the interpolation coordinates (though likely
+depending on the interpolation order or other such parameters), and
+not varying from one output array to another.
+
+The following parameters may be used to query the molecule size:
+\begin{verbatim}
+CCTK_INT const interp_molecule_min_m[N_dims];
+CCTK_INT const interp_molecule_max_m[N_dims];
+\end{verbatim}
+The semantics of these are that if these keys are present (the values
+don't matter), then the interpolator will (re)set the values to give
+the (inclusive) minimum and maximum \verb|m|~molecule coordinates.
+
+The following parameter may be used to query the molecule positions:
+\begin{verbatim}
+CCTK_INT *const interp_molecule_positions[N_dims];
+\end{verbatim}
+The semantics of this is that the caller should set
+\verb|interp_molecule_positions[]| to an array of
+\verb|N_dims| pointers to (caller-supplied) arrays of
+\verb|N_interp_points| \verb|CCTK_INT|s each.  If this key exists,
+then the interpolator will store the molecule positions in the
+pointed-to arrays.
+
+The following parameters may be used to query the
+Jacobian~\eqref{eqn-Jacobian} itself:
+\begin{verbatim}
+CCTK_REAL *const Jacobian_pointer[N_output_arrays];
+const CCTK_INT   Jacobian_offset [N_output_arrays];
+
+/* the next 2 table entries are shared by all Jacobians */
+const CCTK_INT Jacobian_interp_point_stride;
+const CCTK_INT Jacobian_m_strides[N_dims];
+\end{verbatim}
+Then the interpolator would store the Jacobian~\eqref{eqn-Jacobian} in
+\begin{verbatim}
+Jacobian_pointer[out][offset + pt*Jacobian_interp_point_stride
+                             + mi*stride_i + mj*stride_j + mk*stride_k]
+\end{verbatim}
+where
+\begin{verbatim}
+offset = Jacobian_offset[out]
+(stride_i,stride_j,stride_k) = Jacobian_m_strides[]
+\end{verbatim}
+
+By appropriately setting the various stride parameters, this allows
+a fairly wide variety of possible storage layouts for the Jacobian.
+
+An example may help to clarify this:  Suppose we have a 1-D grid
+with 11~grid points, with integer subscripts~0 through~10 inclusive,
+and interpolation coordinates given by \verb|coord_origin = 0.0|
+and \verb|coord_delta = 0.1|.
+Suppose further that we're doing quadratic interpolation, using
+3-point (vertex-centered) molecules, and that the interpolator
+centers the interpolation molecules as close to the interpolation
+point as possible subject to the constraint that the molecules
+never require data from outside the grid.
+Finally, suppose that we query the Jacobian molecule positions for
+the \verb|N_interp_points=14| interpolation points 0.0, 0.04, 0.06,
+0.10, 0.14, 0.16, 0.20, 0.80, 0.84, 0.86, 0.90, 0.94, 0.96, 1.00.
+Then the queries might return
+\begin{verbatim}
+interp_molecule_min_m = -1
+interp_molecule_max_m = +1
+                                         /* interp_x      molecule     */
+interp_molecule_positions[0][ 0] =  1    /*   0.00     [0.0, 0.1, 0.2] */
+interp_molecule_positions[0][ 1] =  1    /*   0.04     [0.0, 0.1, 0.2] */
+interp_molecule_positions[0][ 2] =  1    /*   0.06     [0.0, 0.1, 0.2] */
+interp_molecule_positions[0][ 3] =  1    /*   0.10     [0.0, 0.1, 0.2] */
+interp_molecule_positions[0][ 4] =  1    /*   0.14     [0.0, 0.1, 0.2] */
+interp_molecule_positions[0][ 5] =  2    /*   0.16     [0.1, 0.2, 0.3] */
+interp_molecule_positions[0][ 6] =  2    /*   0.20     [0.1, 0.2, 0.3] */
+interp_molecule_positions[0][ 7] =  8    /*   0.80     [0.7, 0.8, 0.9] */
+interp_molecule_positions[0][ 8] =  8    /*   0.84     [0.7, 0.8, 0.9] */
+interp_molecule_positions[0][ 9] =  9    /*   0.86     [0.8, 0.9, 1.0] */
+interp_molecule_positions[0][10] =  9    /*   0.90     [0.8, 0.9, 1.0] */
+interp_molecule_positions[0][11] =  9    /*   0.94     [0.8, 0.9, 1.0] */
+interp_molecule_positions[0][12] =  9    /*   0.96     [0.8, 0.9, 1.0] */
+interp_molecule_positions[0][13] =  9    /*   1.00     [0.8, 0.9, 1.0] */
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Caching}
+\label{sect-generic-options/caching}
+
+Some interpolators may support special ``caching'' optimizations to
+speed repeated interpolations where some or all of the interpolator
+arguments and/or parameters are the same.  For example, when interpolating
+a tabulated equation of state the number of dimensions, the coordinate
+origin and grid spacing, and the input arrays (the tabulated equation
+of state data), will probably all be the same from one interpolator
+call to another.
+
+If an interpolator supports caching, the following parameters should
+be used to control this:
+
+\begin{verbatim}
+const char cache_type[];        /* set with Util_TableSetString() */
+const char cache_control[];     /* set with Util_TableSetString() */
+CCTK_INT cache_handle;
+\end{verbatim}
+
+There are three basic operations supported:
+\begin{description}
+\item[Create a Cache]
+	To set up caching, call the interpolator with \verb|cache_type|
+	set to describe what arguments and/or parameters will remain
+	the same in future interpolator calls, \verb|cache_control|
+	set to the string \verb|"create"|, and \verb|cache_handle|
+	{\em not\/} in the parameter table.  The interpolator will
+	then do extra (possibly quite time-consuming) work to set
+	up cached information.  The interpolator will delete the
+	key \verb|cache_control|, and return a handle to the cached
+	information in \verb|cache_handle|; this allows multiple caches
+	to be active concurrently.
+\item[Use a Cache]
+	To use a cache (\ie{}~to make an interpolation with the
+	hoped-for speedup), just call the interpolator with
+	\verb|cache_handle| set to the value returned when the cache
+	was created.  Note that you still have to provide all the
+	``will be the same'' interpolator arguments and/or parameters;
+	providing a cache handle is essentially just a promise that
+	these will be the same as in the cache-create interpolator
+	call.  The details of what information is cached, and if/how
+	the ``will be the same'' arguments are still used, are up to
+	the interpolator.
+\item[Destroy a Cache]
+	To destroy a cache (\ie{} free any memory allocated when
+	the cache was created), call the interpolator with
+	\verb|cache_handle| set to the value returned when the cache
+	was created, and  \verb|cache_control| set to the string
+	\verb|"destroy"|.  The interpolator will delete the keys
+	\verb|cache_handle| and \verb|cache_control|, and destroy
+	the cache.
+\end{description}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{The Uniform Cartesian Interpolator}
+
+The uniform Cartesian interpolator is the one which used to live
+in the PUGHInterp thorn.  It was written in C by Thomas Radke in
+early 2001, drawing on earlier Fortran code by Paul Walker.  It
+supports the \verb|CCTK_InterpLocal()| API, and thus doesn't use
+any parameter table at all.  It provides Lagrange polynomial
+interpolation of orders 1, 2, and 3, registered under the operator
+names \verb|"first-order uniform Cartesian"|,
+\verb|"second-order uniform Cartesian"|, and
+\verb|"third-order uniform Cartesian"| respectively.  Each of these
+supports 1, 2, and 3-D arrays.  The code is hard-wired for
+hypercube-shaped interpolation molecules, but it would be easy
+to add additional orders and/or dimensions if desired.
+
+Although the \verb|CCTK_InterpLocal()| API supports both uniform
+and nonuniform grids for the input data, the present implementation
+assumes a uniform grid (and silently gives wrong results for a
+nonuniform grid).
+
+See the Cactus User's Guide ``Full Description of Functions''
+appendix for a full description of the \verb|CCTK_InterpLocal()|
+API, and some examples of how to use it.
+
+%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Implementation}
+
+Internally, this interpolator always does 3-D interpolation, inserting
+zero coordinates as appropriate for lower dimensionalities.  The
+interpolation is done by successive 1-D interpolations along each
+axis.  See the \verb|README| file in the source code directory
+\verb|LocalInterp/src/UniformCartesian/| for further details.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{The Generalized Polynomial Interpolator}
+
+The generalized polynomial interpolator was written in C
+(plus Maple to generate the coefficients) by Jonathan Thornburg in
+winter 2001-2002.  It provides Lagrange polynomial interpolation of
+orders 1-6 for 1-D arrays, and 1-4 for 2- and 3-D
+arrays, all registered under the operator name
+\verb|"generalized polynomial"|.  (Again, it would be easy to add
+additional orders and/or dimensions if desired.)  The code allows
+arbitrarily-shaped interpolation molecules, but at the moment only
+hypercube-shaped molecules are implemented.
+
+This interpolator supports a number of the options described in
+section~\ref{sect-generic-options}:
+\begin{itemize}
+\item	interpolation order
+\item	handling of out-of-range interpolation points
+	(if there are multiple out-of-range points/axes, the
+	one reported will be the first, \ie{}~the out-of-range
+	point with the smallest \verb|pt|, and of that points
+	out-of-range axes, the one with the smallest \verb|axis|)
+	
+\item	non-contiguous input arrays
+\item	derivatives
+	(when taking derivatives with this interpolator,
+	it's most efficient to group all operations on
+	a given input array together in the \verb|operand_indices|
+	and \verb|operation_codes| arrays, as in the example in
+	section~\ref{sect-example-derivatives})
+\end{itemize}
+At present this interpolator does not support computing the Jacobian,
+but we hope to add this in the near future.  This interpolator also
+does not support any caching features; at present (Feb 2002) we have
+no plans to add these.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Smoothing}
+
+The way the generalized polynomial interpolator is implemented it's
+easy to also do Savitzky-Golay smoothing.%%%
+\footnote{%%%
+	 See section 14.8 of Numerical Recipes
+	 2nd edition for a general discussion of
+	 Savitzky-Golay smoothing.
+	 }%%%
+{}  This is best described by way of an example:  Suppose we're doing
+1-D cubic interpolation, \ie{} at each interpolation point we're using
+a cubic polynomial fitted to 4~surrounding data points.  For
+Savitzky-Golay smoothing, we would instead {\em least-squares fit\/}
+a cubic polynomial to some {\em larger\/} number of surrounding data
+points.  This combines interpolation with smoothing, so there's less
+amplification of noise in the input data in the interpolation outputs.
+
+The optional input
+\begin{verbatim}
+const CCTK_INT smoothing;
+\end{verbatim}
+specifies how much (how many points) to enlarge the interpolation
+molecule for this.  The default is 0 (no smoothing).  1 would mean to
+enlarge the molecule by 1 point, \eg{} to use a 5-point molecule instead
+of the usual 4-point one for cubic interpolation.  2 would mean to
+enlarge by 2 points, \eg{} to use a 6-point molecule for cubic
+interpolation.  Etc etc.
+
+Note that in $>1$~dimension, the default hypercube-shaped molecules
+already use more data points than the number of free parameters in
+the interpolating polynomials, \ie{} they already do some Savitzky-Golay
+smoothing.  For example, in 2-D a generic cubic polynomial
+$f(x,y) = \sum_{i+j \leq 3} c_{ij} x^i y^j$ has 10 free parameters
+$c_{ij}$, which we least-squares fit to the 16 data points in the
+$4 \times 4$ molecule.
+
+Savitzky-Golay smoothing is basically free apart from the increase in
+the molecule size, e.g. a \verb|smoothing|=2 cubic interpolation has
+exactly the same cost as any other 6-point--molecule interpolation.
+
+Alas, at the moment only the trivial case \verb|smoothing|=0 is
+implemented, but the framework is all there for more general cases.
+
+%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Implementation}
+
+This interpolator's basic design is to use separate specialized code
+for each combination of
+$(\verb|N_dims|, \verb|molecule_family|, \verb|order|, \verb|smoothing|)$,
+\ie{}~in practice for each distinct choice of interpolation molecule.
+Maple is used to generate all the interpolation coefficients.
+The C preprocessor is then used to generate all the specialized code
+from a single master ``template'' file.  The template code uses
+\verb|#ifdef|s to handle lower dimensions with no extra overhead,
+\eg{}~1-D/2-D interpolation is basically just as efficient as in
+a purpose-built 1-D/2-D interpolator, with no extra overhead imposed
+by the interpolator also supporting higher-dimensional interpolation.
+
+The Maple code which generates the interpolation coefficients is quite
+general-purpose, and can handle an arbitrary dimensionality and molecule
+size/shape.  Generating new coefficients can be rather time-consuming,
+though, \eg{}~the current coefficients for 3-D for orders~1-4 take about
+8~cpu minutes to generate using Maple~7 on a 1.7~GHz~P4.
+
+See the \verb|README| file in the source code directory
+\verb|LocalInterp/src/UniformCartesian/| for further details on the
+implementation.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{A Simple Example of {\tt CCTK\_InterpLocalUniform} Usage}
+
+Here's a simple example of interpolating a \verb|CCTK_REAL| and a
+\verb|CCTK_COMPLEX| $10 \times 20$ 2-D array, at 5 interpolation points,
+using cubic interpolation.
+
+\begin{verbatim}
+#define N_DIMS   2
+#define N_INTERP_POINTS   5
+#define N_INPUT_ARRAYS    2
+#define N_OUTPUT_ARRAYS   2
+
+/* (x,y,z) coordinates of data grid points */
+#define X_ORIGIN   ...
+#define X_DELTA    ...
+#define Y_ORIGIN   ...
+#define Y_DELTA    ...
+const CCTK_REAL origin[N_DIMS] = { X_ORIGIN, Y_ORIGIN };
+const CCTK_REAL delta [N_DIMS] = { X_DELTA,  Y_DELTA  };
+
+/* (x,y,z) coordinates of interpolation points */
+const CCTK_REAL interp_x[N_INTERP_POINTS];
+const CCTK_REAL interp_y[N_INTERP_POINTS];
+const void *const interp_coords[N_DIMS]
+        = { (const void *) interp_x, (const void *) interp_y };
+
+/* input arrays */
+/* ... note Cactus uses Fortran storage ordering, i.e. X is contiguous */
+#define N_X   10
+#define N_Y   20
+const CCTK_REAL    input_real   [N_Y][N_X];
+const CCTK_COMPLEX input_complex[N_Y][N_X];
+const CCTK_INT input_array_dims[N_DIMS] = { N_X, N_Y };
+const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_COMPLEX };
+const void *const input_arrays[N_INPUT_ARRAYS]
+        = { (const void *) input_real, (const void *) input_complex };
+
+/* output arrays */
+CCTK_REAL    output_real   [N_INTERP_POINTS];
+CCTK_COMPLEX output_complex[N_INTERP_POINTS];
+const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_COMPLEX };
+void *const output_arrays[N_OUTPUT_ARRAYS]
+        = { (void *) output_real, (void *) output_complex };
+
+int operator_handle, param_table_handle;
+operator_handle = CCTK_InterpHandle("my interpolation operator");
+if (operator_handle < 0)
+        CCTK_WARN(-1, "can't get interpolation handle!");
+param_table_handle = Util_TableCreateFromString("order=3");
+if (param_table_handle < 0)
+        CCTK_WARN(-1, "can't create parameter table!");
+
+if (CCTK_InterpLocalUniform(N_DIMS,
+                            operator_handle, param_table_handle,
+                            origin, delta,
+                            N_INTERP_POINTS,
+                               CCTK_VARIABLE_REAL,
+                               interp_coords,
+                            N_INPUT_ARRAYS,
+                               input_array_dims,
+                               input_array_type_codes,
+                               input_arrays,
+                            N_OUTPUT_ARRAYS,
+                               output_array_type_codes,
+                               output_arrays) < 0)
+        CCTK_WARN(-1, "error return from interpolator!");
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{An Example of Interplating Derivatives}
+\label{sect-example-derivatives}
+
+Consider the problem described earlier: the computation of all the
+2nd~derivatives of the 3-metric at a set of interpolation points on
+a 2-sphere.  Here's how we could do this in C:
+\begin{verbatim}
+#define N_DIMS        3
+
+/* interpolation points */
+#define N_INTERP_POINTS 1000
+const CCTK_REAL interp_x[N_INTERP_POINTS],
+                interp_y[N_INTERP_POINTS],
+                interp_z[N_INTERP_POINTS];
+const void *const interp_coords[N_DIMS]
+        = { (const void *) interp_x,
+            (const void *) interp_y,
+            (const void *) interp_z };
+
+/* dimensions of the data grid */
+#define NX   30
+#define NY   40
+#define NZ   50
+
+/* input arrays */
+#define N_INPUT_ARRAYS  6
+const CCTK_REAL gxx[NX][NY][NZ], gxy[NX][NY][NZ], gxz[NX][NY][NZ],
+                                 gyy[NX][NY][NZ], gyz[NX][NY][NZ],
+                                                  gzz[NX][NY][NZ];
+
+const CCTK_INT input_array_dims[N_DIMS] = {NX, NY, NZ};
+const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                                    CCTK_VARIABLE_REAL };
+const void *const interp_coords[N_INPUT_ARRAYS]
+        = { (const void *) gxx, (const void *) gxy, (const void *) gxz,
+                                (const void *) gyy, (const void *) gyz,
+                                                    (const void *) gzz };
+
+/* output arrays */
+#define N_OUTPUT_ARRAYS 36
+CCTK_REAL
+ dxx_gxx[N_INTERP_POINTS], dxy_gxx[N_INTERP_POINTS], dxz_gxx[N_INTERP_POINTS],
+                           dyy_gxx[N_INTERP_POINTS], dyz_gxx[N_INTERP_POINTS],
+                                                     dzz_gxx[N_INTERP_POINTS],
+ dxx_gxy[N_INTERP_POINTS], dxy_gxy[N_INTERP_POINTS], dxz_gxy[N_INTERP_POINTS],
+                           dyy_gxy[N_INTERP_POINTS], dyz_gxy[N_INTERP_POINTS],
+                                                     dzz_gxy[N_INTERP_POINTS],
+ dxx_gxz[N_INTERP_POINTS], dxy_gxz[N_INTERP_POINTS], dxz_gxz[N_INTERP_POINTS],
+                           dyy_gxz[N_INTERP_POINTS], dyz_gxz[N_INTERP_POINTS],
+                                                     dzz_gxz[N_INTERP_POINTS],
+ dxx_gyy[N_INTERP_POINTS], dxy_gyy[N_INTERP_POINTS], dxz_gyy[N_INTERP_POINTS],
+                           dyy_gyy[N_INTERP_POINTS], dyz_gyy[N_INTERP_POINTS],
+                                                     dzz_gyy[N_INTERP_POINTS],
+ dxx_gyz[N_INTERP_POINTS], dxy_gyz[N_INTERP_POINTS], dxz_gyz[N_INTERP_POINTS],
+                           dyy_gyz[N_INTERP_POINTS], dyz_gyz[N_INTERP_POINTS],
+                                                     dzz_gyz[N_INTERP_POINTS],
+ dxx_gzz[N_INTERP_POINTS], dxy_gzz[N_INTERP_POINTS], dxz_gzz[N_INTERP_POINTS],
+                           dyy_gzz[N_INTERP_POINTS], dyz_gzz[N_INTERP_POINTS],
+                                                     dzz_gzz[N_INTERP_POINTS];
+const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS]
+   = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                               CCTK_VARIABLE_REAL,
+       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                               CCTK_VARIABLE_REAL,
+       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                               CCTK_VARIABLE_REAL,
+       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                               CCTK_VARIABLE_REAL,
+       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                               CCTK_VARIABLE_REAL,
+       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                               CCTK_VARIABLE_REAL };
+void *const output_arrays[N_OUTPUT_ARRAYS]
+   = { (void *) dxx_gxx, (void *) dxy_gxx, (void *) dxz_gxx,
+                         (void *) dyy_gxx, (void *) dyz_gxx,
+                                           (void *) dzz_gxx,
+       (void *) dxx_gxy, (void *) dxy_gxy, (void *) dxz_gxy,
+                         (void *) dyy_gxy, (void *) dyz_gxy,
+                                           (void *) dzz_gxy,
+       (void *) dxx_gxz, (void *) dxy_gxz, (void *) dxz_gxz,
+                         (void *) dyy_gxz, (void *) dyz_gxz,
+                                           (void *) dzz_gxz,
+       (void *) dxx_gyy, (void *) dxy_gyy, (void *) dxz_gyy,
+                         (void *) dyy_gyy, (void *) dyz_gyy,
+                                           (void *) dzz_gyy,
+       (void *) dxx_gyz, (void *) dxy_gyz, (void *) dxz_gyz,
+                         (void *) dyy_gyz, (void *) dyz_gyz,
+                                           (void *) dzz_gyz,
+       (void *) dxx_gzz, (void *) dxy_gzz, (void *) dxz_gzz,
+                         (void *) dyy_gzz, (void *) dyz_gzz,
+                                           (void *) dzz_gzz };
+
+/* integer codes to specify the derivatives */
+/* (for best efficiency we group all operations on a given input together) */
+const CCTK_INT operand_indices[N_OUTPUT_ARRAYS]
+   = { 0, 0, 0, 0, 0, 0,
+       1, 1, 1, 1, 1, 1,
+       2, 2, 2, 2, 2, 2,
+       3, 3, 3, 3, 3, 3,
+       4, 4, 4, 4, 4, 4,
+       5, 5, 5, 5, 5, 5 };
+#define DERIV(x)   x
+const CCTK_INT operation_codes[N_OUTPUT_ARRAYS]
+   = { DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
+       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
+       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
+       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
+       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
+       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33) };
+
+int operator_handle, param_table_handle;
+operator_handle = CCTK_InterpHandle("my interpolation operator");
+if (operator_handle < 0)
+        CCTK_WARN(-1, "can't get interpolation handle!");
+
+param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT);
+if (param_table_handle < 0)
+        CCTK_WARN(-1, "can't create parameter table!");
+if (Util_TableSetInt(param_table_handle, 3, "order") < 0)
+        CCTK_WARN(-1, "can't set order in parameter table!");
+if (Util_TableSetIntArray(param_table_handle,
+                          N_OUTPUT_ARRAYS, operand_indices,
+                          "operand_indices") < 0)
+        CCTK_WARN(-1, "can't set operand_indices array in parameter table!");
+if (Util_TableSetIntArray(param_table_handle,
+                          N_OUTPUT_ARRAYS, operation_codes,
+                          "operation_codes") < 0)
+        CCTK_WARN(-1, "can't set operation_codes array in parameter table!");
+
+if (CCTK_InterpLocalUniform(N_DIMS,
+                            operator_handle, param_table_handle,
+                            origin, delta,
+                            N_INTERP_POINTS,
+                               CCTK_VARIABLE_REAL,
+                               interp_coords,
+                            N_INPUT_ARRAYS,
+                               input_array_dims,
+                               input_array_type_codes,
+                               input_arrays,
+                            N_OUTPUT_ARRAYS,
+                               output_array_type_codes,
+                               output_arrays) < 0)
+        CCTK_WARN(-1, "error return from interpolator!");
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%
+% Automatically created from the ccl files 
+% Do not worry for now.
+%
+\include{interface}
+\include{param}
+\include{schedule}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Acknowledgments}
+
+Thanks to Ian Hawke for helpful comments on the documentation,
+to Tom Goodale and Thomas Radke for many useful design discussions,
+and to all the Cactus crew for a great infrastructure!
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\end{document}