This file provides basic instructions for use of the Melkes package of
interpolation routines.  It is brief and insufficient, and will be
revised as soon as the author returns from some travel, later in
October, 2001.

Please see first the companion README file for an overview of the
scope of the package.

There are nine core routines, forming in a sense a 3*3 matrix.  Their
names follow the pattern melkes[123]_eval[fgh].

Routines melkes1_eval[fgh] are associated with Melkes interpolation of
parameter p=1, which is just multilinear interpolation.  This case
requires only function values on the vertices of a box.

Routines melkes2_eval[fgh] are for Melkes interpolation of parameter
p=2, which is the cubic case in the sense that the space of
interpolating polynomials contains all monomials that are of total
degree 3 or less.  (Except in the case of dimension one it also
contains some polynomials of total degree greater than 3, but not a
complete set.)  This case requires function values and first
derivatives on the vertices of a box.

Routines melkes3_eval[fgh] are for Melkes interpolation of parameter
p=3, which is the quintic case in the same sense as above.  This case
requires function values and first and second derivatives on the
vertices of a box.

The suffix f, g or h denotes the required calculation.  Routines
melkes[123]_evalf evaluate the interpolant only, at a single point.
Routines melkes[123]_evalg evaluate the interpolant and its gradient
at a single point, and routines melkes[123]_evalh evaluate the
interpolant and its gradient and Hessian at a single point.

The routines melkes[123]_eval[fgh] are general with respect to the
dimensionality of the problem.  In addition, for reasons of
computational efficiency, there are specialized routines for the cases
of dimension in the range 1..4.  These routines have names according
to the pattern melkes[123]_eval[fgh]_d[1234], where the final digit
denotes the dimension.

The routines melkes[123]_eval[fgh] all carry out a certain amount of
error checking, in particular testing that the dimensions of all array
arguments are properly matched.  This error checking is bypassed in
the case of a direct call to the *_d[1234] routines.  For dimensions
greater than 4, only the call to the general routine with no dimension
suffix is available.

It is worth emphasizing at this point that the routines in the Melkes
package deal on each invocation with just a single cell (or box) in
d-dimensional space.  So, if an interpolant is wanted at some point
then the identification of the appropriate cell within the caller's
data structure falls outside the scope of these routines.  The calling
program finds the appropriate cell, and passes to the Melkes routines
only data corresponding to that cell.  It is not even required that
the point x at which the interpolant is evaluated lie inside that
cell; if it lies outside then the interpolating polynomial is still
perfectly well defined and will be evaluated without fuss.

The calling sequence and all arguments for the melkes[123]_eval[fgh]
routines are carefully described in the header sections of the source
code, especially in the *_evalf routine headers.  The reader and
prospective users is, for now, referred to those header sections; they
will not be repeated here.

It will be seen that the melkes1 routines require a real array
argument f0(0:nv-1), which supplies the function values at all
vertices.

The melkes2 routines require arguments f0(0:nv-1), with the same
meaning as in melkes1, and f1(0:nd-1,0:nv-1).  Here, nd is the number
of dimensions (and nv=2**nd is the number of vertices).  The argument
f1 is related quite closely to the gradients at vertices, but it is
not exactly the gradients.  The melkes3 routines require an additional
argument f2(0:nd*(nd+1)/2-1,0:nv-1), which is closely related to the
Hessian but is not the same as the Hessian.

These coefficient arrays f1 and f2 may be computed in routines
melkes2_coef and melkes3_coef, which take as input the proper function
values, gradients and Hessians on the vertices of the box under
consideration.  Routine melkes2_coef takes as input function values
and gradients (f0 and gf) returning the coefficient array f1, and
routine melkes3_coef takes as input function values, gradients and
Hessians (f0, gf and hf) returning the two coefficient arrays f1 and
f2.  The call sequence and arguments to these routines are also
carefully described in the routine headings and that description will
not be repeated here.

The package comes with an extensive set of procedures to facilitate
testing and timing.

Routines melkes[123]_test are testing routines, melkes[123]_time are
timing routines, and melkes[123]_tfun provide an example test function
for the package.  In each of the three cases the test function in
melkes[123]_tfun is of the maximal order for which the interpolation
is mathematically exact, and the test function is devoid of any
obvious symmetries.  The reader is again rferred to the routine header
sections for description of call sequence and arguments.

The package is organized into Fortran-90 modules, of which the
outermost ones are melkes.mod and melkes_test.mod .  Module melkes.mod
combines modules melkes[123].mod, each of which is devoted to one of
the three implemented classes of melkes interpolation.  In turn, each
of the melkes[123].mod modules combines, generally, a 32-bit real and
a 64-bit real version to provide a generic interface.  The testing
module melkes_test.mod likewise combines two precisions.  Specific
modules for one precision would be available as, for example,
melkes2_r8.mod, but it is expected that users will generally not have
a need for those specific interfaces.

There is some small work still to be done on the package content; in
particular I will supply some interfaces to facilitate work with data
on a mesh in the specific cases of diensions two and three.

Bas Braams, Courant Institute, NYU -- October 8, 2001.