Princeton Splines -- D. McCune Feb. 1999, updated Feb. 2004.
(software contributions also by W. Houlberg and A. Pletzer).
This document describes PSPLINE -- a library of Spline and Hermite
cubic interpolation routines for 1d, 2d, and 3d datasets on
rectilinear grids.
Abstract --
PSPLINE -- a collection of Spline and Hermite interpolation tools for
1D, 2D, and 3D datasets on rectilinear grids. Spline routines give
full control over boundary conditions -- the user may specify "periodic",
"not a knot", 1st derivative match, or 2nd derivative match boundary
conditions on either end of each grid dimension. Hermite routines take
as input the function value and derivatives at each grid point, giving
back a representation of the function between grid points. The splines
are continuously twice differentiable in all directions across all
grid cell boundaries and over the entire grid domain. Hermite functions
are continuously once differentiable in all directions over the entire
grid domain. For a representation of dimensionality N, an N-dimensional
Hermite or Spline function requires 2**N*(nx1*nx2*...*NxN) memory
words; there is also a somewhat faster but much more memory intensive
Spline representation requiring 4**N*(nx1*nx2*...*nxN) memory words.
Speed of evaluation of Hermite and Spline functions are similar--
both are constructed of cubic polynomial segments or "patches", which
are often much faster to evaluate than representations using Fourier
series or otherwise involving transcendental functions.
Both spline and hermite methods offer highly efficient options for
differentiable representations of higher dimensional objects-- see the
case study.
A defect in a 1d spline fitting routine caused some spline fits using the
higher order boundary conditions
d2f/dx2 .ne. 0
d2f/dx2 from 2nd divided difference (edge parabolic fit)
d3f/dx3 from 3rd divided difference (edge cubic fit)
To be C0 instead of C2. This has been repaired.
The compact bicubic and tricubic splines are now constructed directly,
rather than by conversion from the non-compact representations. In the
case of the tricubics, this saves a factor of nearly 8 on memory usage
and a factor of about 2-5 on cpu time for construction of the spline
coefficients, depending on problem size, workstation cache size, etc.
I would like to acknowledge Kenneth Esler, a graduate student in
condensed matter physics at University of Illinois -- Urbana Champaign,
who downloaded pspline from the NTCC modules library.
In the course of his use of tricubic splines in his condensed matter
physics application, he encountered performance problems due to the
memory-intensive "non-compact" construction method used in the
previous version of the code. He sent a C++ code sample suggesting
a more compact construction method. I found his suggestion to be
usable; a fortran-90 implementation for compact construction of
bicubic and tricubic splines was developed and is now included in
the library. The fortran-90 interface (ezspline) uses this new method.
In the fortran-77 interface, use:
subroutine mkbicub -- construct compact bicubic C2 spline
subroutine mktricub -- construct compact tricubic C2 spline
These routines were completed on Feb. 26 2004 -- by Douglas McCune.
The advice of PPPL physicist Dr. Leonid Zakharov is much appreciated,
for pointing out the existence of the 2**N "compact" spline
representation. Dr. Zakharov and others also pointed out the
desirability of including cross derivative terms in the higher
dimensional Hermite formulations.
I also want to express thanks to ORNL physicist Dr. Wayne Houlberg
for reviewing the PSPLINE software and documentation, and for
writing "v_spline" which re-implemented the 1d spline tridiagonal
solver in a cleaner, more efficient form, with additional boundary
condition options.
Finally, I want to express appreciation to PPPL Computational
Physicist Dr. Alex Pletzer for his contribution of the EZspline
software interface to pspline.
PPPL Computational Plasma Physics Group, Feb. 5 1999 -- D. McCune
Updated: June 10, 1999
A 3d Spline Representation of B field greatly speeds stellarator fast
ion orbit calculations --
CPPG has developed two sets of highly general tricubic spline
interpolation tools for 3d objects. Periodic boundary condition
options make the 3d splines especially well suited to representation
of objects such as B(psi,theta,zeta) which are functions of poloidal
and toroidal angle parameters which arise in magnetic flux coordinate
systems. The spline tools have now been installed in R. White's
ORBIT code [1,2].
Two types of tricubic splines are supported:
* "explicit" spline-- the fastest, but very memory intensive:
64*Npsi*Ntheta*Nzeta memory words for the spline coefficients
* "compact" spline-- not quite as fast but eight times less
memory consumption: 8*Npsi*Ntheta*Nzeta memory words.
Evaluation of the 3d spline, which is continuously twice differentiable
in all coordinates, consists (for the explicit spline) in computing a
sum of 21 cubic polynomials-- regardless of the complexity of the object
represented (however, more complex objects require finer meshes and much
more memory). The fixed cost of evaluation distinguishes the spline
method from the Fourier harmonics summation representation typically
used in the past, which requires additional harmonics to represent
more complicated objects or fields, making each evaluation
proportionately more expensive.
In a trial case-- a stellarator equilibrium B field representation
requiring 304 Fourier harmonics to describe the poloidal and toroidal
variation-- replacement of the Fourier representation with a 3d spline
representation resulted in evaluation speed 80x faster on a DEC alpha.
By controlling the fineness of the mesh, the spline can reproduce the
Fourier representation to any desired accuracy. On the other hand,
the memory cost of the spline is significantly higher than the moments
representation, which needed just 4 x Nmoments x Npsi memory words.
The 3d spline has been vectorized and installed and tested in Roscoe
White's 3d Stellarator fast ion orbit code ORBIT3D on Killeen (a vector
Cray supercomputer at NERSC). For a reference problem (1000 particles
and 100 transits), the following table shows the speedups obtained,
relative to the older harmonics based code ORBITMN:
# of harmonics Speedup factor (ORBIT3D vs. ORBITMN)
======================================================
30 4 times faster
100 10 " "
300 50 " "
A typical 3d stellator equilibrium would require on the order of 100
harmonics (but 300 harmonic cases do arise). It is expected that
ORBIT3D, as a replacement to ORBITMN, will be heavily used to study
fast ion confinement in stellarator B field equilibria.
The 3d spline may have application in other codes as well -- e.g.
particle simulations. Because the data structure has locality in
three dimensions, it may be well suited to a 3d domain decomposition
for parallelization.
[June 10, 1999 addendum] ORBIT3D now allows switching between use
of the "explicit" and the "compact" spline method. Although lookups
using the explicit spline were 1.3 times faster than compact splines
on a standalone test application built on a DEC Alpha workstation,
inside ORBIT3D, using hand inlined vectorized interpolant evaluation
code for both methods, the net difference was only about 10%:
ORBIT3D on Cray-J90, 5000 ptcls, 100 transits
kspline=1 1227 sec (explicit spline)
kspline=2 1374 sec (compact spline) [3].
Probably for most 3d applications the compact spline's factor of 8
savings in memory utilization would outweigh the modest loss in
computational performance.
[1] R. B. White and M. S. Chance, Phys Fluids 27, 2455 (1984)
[2] R. B. White, Phys Fluids B 2, 845 (1990)
[3] R. B. White, priv. communication, June 10, 1999.
Cubic Splines are piecewise cubic functions that pass through a
set of given data points. Cubic splines are continuous and
continuously twice differentiable.
1D Cubic Splines have long been popular interpolation tools. Use of
splines always involves two steps: (1) setup of spline coefficients,
and (2) evaluation of spline at desired target points.
The setup of 1D spline coefficients involves a highly stable numerical
procedure-- solution of a diagonally dominant tri-diagonal matrix [1,2].
The following inputs are required:
the grid: x(i), i=1 to Nx x(i+1).gt.x(i) for all i.
the data: f(i), i=1 to Nx
a boundary condition at x(1)
a boundary condition at x(Nx).
This data uniquely specifies a cubic spline function, which is
output from the spline setup routine as a set of spline coefficients.
There are two choices for representation of the spline coefficients--
an "explicit" representation requiring 4*Nx memory words, and a Hermite-
like "compact" representation which stores just the data and 2nd
derivative at the grid points, requiring 2*Nx memory words.
The "explicit" representation is somewhat faster to use and is probably
the best choice for 1d splines. But, for higher dimensional splines,
these two representations require, respectively (where N is the
dimensionality or rank of the object):
explicit -- 4**N*(Nx1*Nx2*...*NxN)
compact -- 2**N*(Nx1*Nx2*...*NxN)
Thus for N=3, the explicit representation requires 8 times the memory
of the compact representation; this memory cost savings is likely to
outweigh the modestly higher computational cost of using the compact
representation scheme.
In 1d, the setup of a spline using the explicit representation yields
an array C(4,Nx) such that the spline function is evaluated by the
formulae
piecewise cubic spline function value:
s(x)=C(1,i)+(x-x(i))*C(2,i)+(x-x(i))**2*C(3,i)+(x-x(i))**3*C(4,i)
for x(i).le.x.le.x(i+1)
piecewise quadratic spline function 1st derivative:
s'(x)=C(2,i)+2*(x-x(i))*C(3,i)+3*(x-x(i))**2*C(4,i)
for x(i).le.x.le.x(i+1)
piecewise linear spline function 2nd derivative:
s''(x)=2*C(3,i)+6*(x-x(i))*C(4,i)
for x(i).le.x.le.x(i+1)
and satisfying the boundary conditions and the requirements:
s(x(i))=f(i), i = 1 to N
s'(x) is continuous accross all interior x(i)
s''(x) is continuous accross all interior x(i)
The compact representation requires an array D(2,Nx) where D(1,j)
gives the function value at grid point x(j), and D(2,j) gives the
interpolant's 2nd derivative at x(j). Then, evaluation takes the
following form.
If x(i).le.x.le.x(i+1), let
h = (x(i+1)-x(i))
p = (x-x(i))/h
q = 1-p
s(x)= q*D(1,i) + p*D(1,i+1) +
(h**2/6)*{[q**3-q]*D(2,i) + [p**3-p]*D(2,i+1)}
s'(x) = (D(1,i+1) - D(1,i))/h +
(h/6)*{[3*p**2-1]*D(2,i+1) - [3*q**2-1]*D(2,i)}
s''(x) = q*D(2,i) + p*D(2,i+1)
where (by the choice of D) the given data points are interpolated and
the specified boundary conditions are satisfied.
Typical boundary conditions given are s'(x(1)) or s''(x(Nx)). Many
spline routines also allow the boundary condition to be inferred from
the data in the vicinity of the boundary, e.g. match a numerical
third derivative from the first four data points as in "spline.for" [1]
or a "not a knot" boundary condition forcing the third derivative to be
continuous at x(2) and x(N-1) as in the default usage of "cubspl.for" [2].
(Note that the third derivative of a cubic spline is a step function
and generally not continuous -- imposition of continuity near the
boundary constitutes the "not a knot" boundary condition).
Another useful option is the "periodic" boundary condition
s'(x(1))=s'(x(Nx))
s''(x(1))=s''(x(Nx)).
(note s(x(1))=s(x(Nx)), i.e. f(1)=f(Nx), is not required).
This requires a simple modification to the spline coefficients
algorithm (typically left as an exercise in the texts) as the
coefficients solution matrix is no longer strictly tridiagonal
(it has non-zero off diagonal corner elements). The PSPLINE
library originally took its name from "pspline.for" which contains
an algorithm for producing a periodic spline.
v_spline.for gathers together all the 1d spline coefficient calculations
into a single routine.
references --
[1] Forsythe, Malcolm, and Moler, "Computer Methods for Mathematical
Computations", Prentice-Hall, 1977. (see chapter on interpolation).
[2] C. De Boor, "A Practical Guide to Splines".
[3] W. Press, S. Teukolsky, W. Vetterling, B. Flannery, "Numerical
Recipies in Fortran 77 (2nd edition)", pp 107-110 -- discusses the
compact representation.
Also useful:
[4] Engeln-Muellges, Uhlig, Numerical Algorithms with Fortran, Springer,
1996, p.251
It turns out that 1D splines can be extended to define interpolation
of data given on a 2d rectilinear grid:
f(i,j), i=1 to Nx, j=1 to Ny,
vs. x(i), i=1 to Nx, x(i).lt.x(i+1), and
vs. y(j), j=1 to Ny, y(j).lt.y(j+1).
As in the 1d case, boundary condition data is required: typically,
df/dx or d2f/dx2 vs. y(j) at x(1) and x(Nx), and, df/dy or d2f/dy2
vs. x(i) at y(1) and Y(Ny). As in the 1d case, sometimes the data
itself can be used to infer a reasonable boundary condition, obviating
the requirement for its explicit specification. Or, perhaps, a periodic
boundary condition can be used. In the case of a bicubic spline
representation of a function f(r,theta) on a polar coordinate system,
the periodic boundary condition would be natural, and the preferred
choice.
For a given choice of Nx*Ny data points f(i,j) and 2*Nx + 2*Ny
boundary conditions, it turns out there is exactly one bicubic
spline function which is continuously twice differentiable over
the entire domain x(1) .le. x .le. x(Nx) & y(1) .le. y .le. y(Ny).
Existence / Uniqueness proofs and a construction method will be
sketched below.
In the PSPLINE library, a choice between two simple representations
is available for bicubic splines. In the "explicit" representation,
for each data point (x(i),y(j)), a 4x4 array of 16 coefficients Ckl
are retained.
The bicubic spline function s(x,y) is then evaluated as follows:
for x(i).le.x.le.x(i+1), let dx=(x-x(i))
for y(j).le.y.le.y(j+1), let dy=(y-y(j))
let C11...C44 denote the 16 coeffcients associated with (i,j).
then:
s(x,y)= C11 + dx*C21 + dx**2*C31 + dx**3*C41
+dy*(C12 + dx*C22 + dx**2*C32 + dx**3*C42)
+dy**2*(C13 + dx*C23 + dx**2*C33 + dx**3*C43)
+dy**3*(C14 + dx*C24 + dx**2*C34 + dx**3*C44)
This can be thought of as the evaluation of 5 cubic pieces, each
of which can be cast in the form A+dx*(B+dx*(C+dx*D)), i.e. three
A*x+B operations: so then a total of 15 A*x+B operations for
a single evaluation of s(x,y)-- thus, a very computationally
efficient representation.
The formulae for ds/dx, ds/dy, d2s/dx2, d2s/dy2, d2s/dxdy are
all derived from the formula for s(x,y) in the obvious way.
The explicit representation requires 16*Nx*Ny memory words. The
compact representation requires only 4*Nx*Ny memory words, and
is only slightly slower computationally. In the compact
representation one uses
D(1,i,j) = the function value at (x(i),y(j))
D(2,i,j) = the 2nd derivative d2s/dx2 at (x(i),y(j))
D(3,i,j) = the 2nd derivative d2s/dy2 at (x(i),y(j))
D(4,i,j) = the 4th derivative d4s/dx2dy2 at (x(i),y(j))
where the spline derivative coefficients are carefully chosen
to achieve the required level of continuity, differentiability,
and satisfaction of boundary conditions.
Evaluation of the spline using the compact representation is
as follows:
for x(i).le.x.le.x(i+1) and y(j).le.y.le.y(j+1), let
hx = x(i+1)-x(i) hy = y(j+1)-y(j)
px = (x-x(i))/hx py = (y-y(j))/hy
qx = 1-px qy = 1-py
rx = (px**3-px) ry = (py**3-py)
sx = (qx**3-qx) sy = (qy**3-qy)
s(x,y) = qx*(qy*D(1,i,j)+py*D(1,i,j+1)) +
px*(qy*D(1,i+1,j)+py*D(1,i+1,j+1)) +
(hx**2/6)*{sx*(qy*D(2,i,j)+py*D(2,i,j+1)) +
rx*(qy*D(2,i+1,j)+py*D(2,i+1,j+1))} +
(hy**2/6)*{qx*(sy*D(3,i,j)+ry*D(3,i,j+1)) +
px*(sy*D(3,i+1,j)+ry*D(3,i+1,j+1))} +
(hx**2*hy**2/36)*
{sx*(sy*D(4,i,j)+ry*D(4,i,j+1)) +
rx*(sy*D(4,i+1,j)+ry*D(4,i+1,j+1))}
Of course the key is to evaluate the set of coefficients {Ckl}
or {D} such that the continuously twice differentiable property is
achieved. The fact that this can be done, and how it is done,
is considered in the subtopics. The subtopics focus on obtaining
{Ckl}; the derivation of {D} from {Ckl} is straight-forward.
A proof is sketched here for the following result:
given: an x grid {x(1),x(2),...,x(Nx)}, x(i).lt.x(i+1),
a y grid {y(1),y(2),...,y(Ny)}, y(j).lt.y(j+1)
a data array f(i,j) specifying the desired
interpolation values,
for each y grid point j:
a boundary condition bcx1(j) at x(1) and
bcx2(j) at x(Nx),
and for each x grid point i:
a boundary contition bcy1(i) at y(1) and
bcy2(i) y(Ny),
then: there exists one and only one bicubic spline
function s(x,y) satisifying
s(x(i),y(j)) = f(i,j)
AND all the boundary conditions
AND s(x,y) continuous in [x(1),x(Nx)]x[y(1),y(Nx)]
AND ds/dx continuous in [x(1),x(Nx)]x[y(1),y(Nx)]
AND ds/dy continuous in [x(1),x(Nx)]x[y(1),y(Nx)]
AND d2s/dx2 continuous in [x(1),x(Nx)]x[y(1),y(Nx)]
AND d2s/dy2 continuous in [x(1),x(Nx)]x[y(1),y(Nx)]
AND d2s/dxdy continuous in [x(1),x(Nx)]x[y(1),y(Nx)].
Theoretical results on splines can be had by persuing methods of
classical analysis -- such as, construction of set of independent
"basis" splines from which any spline can be built. The "basis"
splines discussed here are not efficient for computation, but are
useful as a theoretical tool.
It has been noted that boundary conditions for splines can be
expressed in a variety of ways. For simplicity of this discussion
(and without loss of generality) it will be assumed that the
boundary conditions at x(1) and x(Nx) specify ds/dx at each y(j),
and that the boundary conditions at y(1) and y(Ny) specify ds/dy
at each x(i).
Then, the basis set can be formed as follows. First define the
the 1d "x basis" splines (for i=1 to Nx) s[i](x) on
{x(1),x(2),...,x(Nx)} satisfying
s[i](x(j)) = delta(i,j) (i.e. 1 if i=j, 0 otherwise)
and
ds[i]/dx = 0 at x(1)
and
ds[i]/dx = 0 at x(Nx).
Because of the properties of 1d splines, s[i] exists and is
unique. Then define also s[0](x) satisfying
s[0](x(j)) = 0 for all j
and
ds[0]/dx = 1 at x(1)
and
ds[0]/dx = 0 at x(Nx).
Similarly, define s[Nx+1](x) satisfying
s[Nx+1](x(j)) = 0 for all j
and
ds[Nx+1]/dx = 0 at x(1)
and
ds[Nx+1]/dx = 1 at x(Nx+1).
Again the properties of 1d splines show that s[0](x) and s[Nx+1](x)
exist and are unique.
Similarly, define a set of Ny+2 1d "y basis" splines s[j](y) on
{y(1),y(2),...,y(Ny)} satisfying the precisely analogous conditions
on the values of s[j](y(k)) and ds[j]/dy at y(1) and y(Ny). Again,
all of the 1d "y basis" splines exist and are unique.
From the 1d basis splines construct 2d basis splines s[i,j](x,y) by
the simple formula
s[i,j](x,y) = s[i](x) * s[j](y).
It can readily be seen that s[i,j](x(k),y(l)) = delta(i,k)*delta(j,l)
(for the interior), and ds/dx[0,j](x(1)) = delta(j,l) for the x(1)
boundary, and so on for the x(Nx), y(1) and y(Ny) boundaries.
Because s[i](x) and s[j](y) are continuously twice differentiable
in x and y respectively, it is clear that s[i,j] is continuously
twice differentiable on [x(1),x(Nx)]x[y(1),y(Nx)].
Therefore, the sum
s(x,y) = [sum i=1 to Nx, j=1 to Ny] f(i,j)*s[i,j](x,y)
+ [sum j=1 to Ny] bcx1(j)*s[0,j](x,y)
+ [sum j=1 to Ny] bcx2(j)*s[Nx+1,j](x,y)
+ [sum i=1 to Nx] bcy1(i)*s[i,0](x,y)
+ [sum i=1 to Nx] bcy2(i)*s[i,Ny+1](x,y)
is also continuously twice differentiable on
[x(1),x(Nx)]x[y(1),y(Nx)]. But by the definition of the s[i,j] basis
splines, s(x,y) also satisfies all of the required conditions
s(x(i),y(j)) = f(i,j) (i=1 to Nx, j=1 to Ny)
ds/dx(x(1),y(j)) = bcx1(j)
ds/dx(x(Nx),y(j)) = bcx2(j)
ds/dy(x(i),y(1)) = bcy2(i)
ds/dy(x(i),y(Ny)) = bcy2(i)
This proves existence.
Uniqueness is shown in the usual way, using a proof by contradiction.
Assume there exists another bicubic spline s2(x,y) that satisfies all
the conditions. Then the difference
r(x,y) = s(x,y) - s2(x,y)
must satisfy r(x(i),y(j)) = 0 for all i,j, and also dr/dx = 0 at
x(1) and x(Nx), and dr/dy = 0 at y(1) and y(Ny).
The general bicubic spline r(x,y) can be expressed in the form
r(x,y)= C11 + dx*C21 + dx**2*C31 + dx**3*C41
+dy*(C12 + dx*C22 + dx**2*C32 + dx**3*C42)
+dy**2*(C13 + dx*C23 + dx**2*C33 + dx**3*C43)
+dy**3*(C14 + dx*C24 + dx**2*C34 + dx**3*C44)
where each Ckl is associated with a subregion
[x(k),x(k+1)]x[x(l),x(l+1)] of the gridded space.
Consideration of the restriction of r(x,y) to x(i) or y(j)
quickly leads to the conclusion that C11=C21=C31=C41=C12=C13=C14=0
(because a 1d spline with zero values and zero slope boundary
conditions is itself identically zero).
This leaves in each grid square [x(k),x(k+1)]x[x(l),x(l+1)] a set
of 9 coefficients C22,C32,C42,C23,C33,C43,C24,C34,C44 not yet
demonstrated to be zero. But, consideration of continuity and
differentiability requirements across grid square boundaries
quickly leads to an overdetermined system of linear equations
involving these coefficients, the only solution of which is
that all these coefficients, too, are zero.
Thus r(x,y)=0 everywhere, and s(x,y)=s2(x,y) -- contradicting
the assumption of non-uniqueness.
Thus uniqueness is proven.
Given that there exists a unique continuous twice differentiable
bicubic spline interpolating given data and satisfying given
boundary conditions, how is it to be constructed efficiently?
The bicubic spline existence proof showed a method of construction,
but, it is not practical. Each spline basis function is an object
with size of proportional to (Nx+Ny), and there are (Nx+2)*(Ny+2)
such 2d basis splines to be considered. Any algorithm using basis
splines is at best going to have computational performance scaling
of Nx*Ny*(Nx+Ny). An construction algorithm that scales as Nx*Ny
should be sought first. It can be found.
Before describing this in detail, a digression on boundary condition
"sets" is necessary.
The set of boundary conditions on the spline at x(1) and x(Nx) shall
be denoted as the X-set, and the set at y(1) and y(Ny) shall be denoted
the Y-set.
Boundary condition sets can be divided into two categories:
"homogeneous" and "inhomogeneous". A homogeneous boundary condition
set A has the property that if s1(x,y) satisfies A, then (for any real
number c) c*S1 also satisfies A, and, if both s1(x,y) and s2(x,y)
satisfy A then so also does s1+s2.
The following are examples of homogeneous boundary condition X-sets:
* periodic boundary condition:
ds/dx(x(1),y)=ds/dx(x(Nx),y), d2s/dx2(x(1),y)=d2s/dx2(x(Nx),y)
* any set constructed by taking one condition from column A
and one from column B:
A B
========================== =============================
d3s/dx3(x(2),y) continuous d3s/dx3(x(Nx-1),y) continuous[*]
ds/dx(x(1),y)=0 ds/dx(x(Nx),y)=0
d2s/dx2(x(1),y)=0 d2s/dx2(x(Nx),y)=0
d[k]s/dx[n] = k'th divided d[k]s/dx[k] = k'th divided
difference at x(1) [**] difference at x(n) [**]
[*] these are the "not a knot" boundary conditions, which are
homogeneous.
[**] e.g. for the "first" divided difference at x(1):
ds/dx(x(1),y(j)) = [f(2,j)-f(1,j)]/(x(2)-x(1))
2nd and 3rd divided difference expressions for d2s/dx2 and
d3s/dx3 are also supported (software contributed by Wayne
Houlberg).
* inhomogeneous boundary conditions all involve explicitly setting
a derivative or 2nd derivative to a fixed non-zero value at a
boundary point.
Homogeneous boundary conditions are "nice" because they facilitate
the construction of bicubic splines. To see this, consider a
particular bicubic spline construction problem (non-compact
representation). We are seeking a continuous twice differentiable
spline of the form
s(x,y)= C11 + dx*C21 + dx**2*C31 + dx**3*C41
+dy*(C12 + dx*C22 + dx**2*C32 + dx**3*C42)
+dy**2*(C13 + dx*C23 + dx**2*C33 + dx**3*C43)
+dy**3*(C14 + dx*C24 + dx**2*C34 + dx**3*C44)
where the 16 coefficients {Ckl} are defined separately in each
grid square [x(i),x(i+1)]x[y(j),y(j+1)].
and we are given:
* an x grid {x(1),x(2),...,x(Nx)}, x(i).lt.x(i+1),
* a y grid {y(1),y(2),...,y(Ny)}, y(j).lt.y(j+1)
* a data array f(i,j) specifying the desired
interpolation values,
* The X-set of boundary conditions:
for each y grid point j:
a boundary condition bcx1(j) at x(1) and
bcx2(j) at x(Nx),
* The Y-set of boundary conditions:
for each x grid point i:
a boundary contition bcy1(i) at y(1) and
bcy2(i) y(Ny).
which we can denote as problem S.
If we assume that the Y-set is homogeneous, we can construct the
unique twice differentiable bicubic spline by the following
algorithm:
Algorithm A:
(a) at each y(j), compute the 1d spline through the data points
f(x(1),y(j)), f(x(2),y(j)), ... ,f(x(Nx),y(j)) and satisfying
the (not necessarily linear) boundary conditions at (x(1),y(j))
and (x(Nx),y(j)). In each grid square this yields the
coefficients {C11, C21, C31, C41}. Note that for grid cell
(i,j), C11 is just f(i,j).
(b) using the HOMOGENEOUS Y-set boundary condition at each x(i),
find the 1d splines s1(y), s2(y), s3(y) and s4(y) satisfying:
s1(x(i),y(j))=C11(i,j)
s2(x(i),y(j))=C21(i,j)
s3(x(i),y(j))=C31(i,j)
s4(x(i),y(j))=C41(i,j)
which yields four spline coefficients for each 1d spline.
These 4 x 4 coefficients are precisely the {Ckl} sought.
Because the time for computation of an N point 1d spline, (which
is just the solution of a tridiagonal system) is proportional to c*N,
the above algorithm scales 5*c*Nx*Ny -- i.e. proportional to Nx*Ny,
as desired.
Why does this method work? It can be seen by thinking (again) in
terms of basis splines. Consider a simplified problem P[j0] for which
(1) the X boundary conditions and function values are zero for
all y(j) except when j=j0, and
(2) The Y-set of boundary conditions is linear.
Let s[j0](x) be the 1d spline at y(j0) that satisfies the given data
and boundary conditions along along y(j0). Then, there is a 1d basis
spline b(y) which satisfies the lineary Y-set boundary conditions, and,
b(y(j)) = 1 if j=j0
= 0 otherwise.
The solution to P[j0] is just
s(x,y) = s[j0](x)*b(y),
and this solution IS CONSTRUCTED by algorithm A, because all of the
1d y spline problems presented yield solutions proportional to b(y),
and because of the Y-set boundary condition linearity allows these
solutions to be combined with preservation of the boundary condition.
The Y-set linearity also allows the full problem S to be considered
as a sum of P-type problems. It is easily shown that
sum[j0=1 to Ny] A(P[j0]) = A(sum[j0=1 to Ny] P(j0)) = A(S).
=================
Thus, the general problem of a construction algorithm with Nx*Ny
performance scaling is solved, when granted a HOMOGENEOUS Y-set of
boundary conditions.
If we have a problem S' where the Y-set is inhomogeneous, we proceed as
follows (algorithm B):
Let A' be the "complementary" solution algorithm: it solves the
problem with a linear X-set by splining first in y, then splining
the resulting y 1d spline coefficients in x.
Algorithm B involves a combination of algorithms A and A'. To solve
problem S using B:
(a) use algorithm A to solve the related problem S* for a bicubic
spline that satisfies the X-set boundary conditions, and
"not-a-knot" homogeneous boundary condition at y(1) and y(Ny).
(b) use algorithm A' to solve the residual problem R=S'-S*. Note
that the data values and X-set boundary conditions for R are
all zero. Problem R's Y-set boundary conditions are the
difference between the S* boundary condition obtained, and the
desired (inhomogeneous) boundary condition of the original problem.
(c) The sum of the solutions A(S*) and A'(R) is the solution to
the original problem, S'.
The computational cost of this algorithm-- approximately 1.5x the
cost of algorithm A alone-- and scales like (Nx*Ny), and is
acceptable.
The spline coefficients in the non-compact representation are
closely related to the coefficients of the compact representation,
i.e. for any grid cell in a bicubic spline,
C(1,1) = f
2*C(3,1) = fxx
2*C(1,3) = fyy
4*C(3,3) = fxxyy
Consideration of these relations allows a compact bicubic (or
tricubic) spline construction method to be written in a manner
analogous to the non-compact method (but considerably more
efficient: C(2,2), etc. are not re-splined). The issue of
inhomogeneous boundary conditions can be handled in the same
way as it is for the non-compact construction method.
NAG comparison notes -- D. McCune 21 Feb. 1999 --
Summary: NAG bicubic B-splines are about 4x slower to evaluate,
but, 16x more compact in memory consumption then the "explicit"
representation, and 4x more compact than the "compact"
representation used in PSPLINE.
Tests comparing PSPLINE r8bcspeval and NAG e02def indicate that
r8bcspeval is 4x faster. The B-spline representation created in NAG
e01daf and evaluated in NAG e02def is far more compact -- consuming
sixteen times less memory. But, there are no options for controlling
boundary conditions in the nag bicubic b- spline routine. Also, there
are no 3d routines.
The comments about the extensibility of 1d splines to 2d applies equally
well to 3d rectilinear grids. So, too, do the comments about the method
of construction, the differentiability properties of the result, and
the choice between the (faster) "explicit" representation requiring
64*Nx*Ny*Nz memory words, and the "compact" representation requiring
8*Nx*Ny*Nz memory words.
The input data for setup of a 3d spline are
f(i,j,k) at x(i),i=1 to Nx, y(j), j=1 to Ny, z(k), i=1 to Nz.
with x(i+1).gt.x(i),
y(j+1).gt.y(j), and
z(k+1).gt.z(k) required.
General boundary conditions are 1st or second derivative information
across a grid boundary plane, e.g. df/dx at (x(1),y(j),z(k)), an Ny x Nz
array of data.
Periodic and "not a knot" boundary conditions are again available. For
data given on a toroidal coordinate system, f(r,theta,phi), periodic
boundary conditions in the theta and phi directions would be the
preferred choice.
The storage requirement for the explicit 3d spline representation
data set is 64*Nx*Ny*Nz. The cost of evaluation is 21 computations
of the form ((A*x+B)*x+C)*x+D (63 A*x+B type operations). The storage
requirement for the compact representation is 8*Nx*Ny*Nz, and timing
tests (on a DEC alpha workstation) show the cost of evaluation to be
approximately 1.4 times greater than for the explicit representation.
The PSPLINE library offers three routines, with different boundary
condition option sets, for the setting up of tricubic spline data sets.
In addition, there is a scalar bicubic spline evaluation routine and
a vectorized spline evaluation routine.
Remarkably, standard mathematics software libraries (NAG, netlib) do not
seem to include tricubic splines. Yet this representation can be highly
efficient in practice -- the author has used a tricubic spline to replace
an optimized hybrid Fourier series & 1d spline representation of an
actual physics model object, with a speedup in evaluation of a factor
of nearly 100 on a DEC alpha workstation (although a substantial memory
cost is incurred). See the "Case Study" (top level subtopic).
Certainly 4d and even 5d multi-cubic splines could be constructed
using the compact representation. However, the storage and evaluation
costs (algebra) are exponential in the dimensionality.
Higher-dimensional splines (beyond 3d) are not currently available in
the library.
Hermite cubic interpolation is sometimes used as an alternative to
splines. The Hermite method avoids the global tridiagonal matrix
inversion for coefficients evaluation required by splines, and the
local behaviour of the Hermite interpolant is controlled strictly
by local data -- sometimes an advantage. However, the Hermite
interpolant is continuously differentiable only once.
Hermite interpolation requires function values and derivatives
to be provided by the user at all grid points. The following
inputs are required:
the grid: x(i), i=1 to N x(i+1).gt.x(i) for all i.
the data: f(0,i), i=1 to N f(0,i)=f(x(i))
the derivative: f(1,i), i=1 to N, f(1,i) = [df/dx](x(i)).
In 1d, the construction of the Hermite interpolant is extremely simple.
In each interval [x(i),x(i+1)], there is only one cubic polynomial
that matches the data values f(0,i), f(0,i+1) and the derivative values
f(1,i), f(1,i+1) at the end points. This cubic polynomial is the
Hermite cubic for that interval. The value and derivative of the
Hermite cubic is continuous with the Hermite cubics in the neighbouring
intervals simply by virtue of the sharing the data value and derivative
information {f(0,i), f(1,i)} and {f(0,i+1),f(1,i+1)}.
In practice, it is not necessary to solve for the coefficients of the
cubic piece in each interval. Rather, it can be constructed on the
fly from the Hermite basis functions, which, on the unit interval [0,1]
are given by:
alpha(x)=x*x*(-2*x + 3)
alpha(0)=0 alpha(1)=1 alpha'(0)=0 alpha'(1)=0
alphabar(x)=1-alpha(x)
alphabar(0)=1 alphabar(1)=0 alphabar'(0)=0 alphabar'(1)=0
beta(x)=x*x*(x-1)
beta(0)=0 beta(1)=0 beta'(0)=0 beta'(1)=1
betabar(x)=x*(x*x-2*x+1)
betabar(0)=0 betabar(1)=0 betabar'(0)=1 betabar'(1)=0
Speed of evaluation of values and derivatives on a general grid requires
access to the information hx(i)=(x(i+1)-x(i)) and hxi(i)=1/hx(i). Given
this information, Hermite evaluation speeds are similar to spline
evaluation speeds.
The memory consumption of Hermite and Spline single and multivariate
representations are equivalent, as per the following table:
1d 2d 3d M-dimensional
Hermite 2*N 4*Nx*Ny 8*Nx*Ny*Nz 2**M*(Nx1*Nx2*...*NxM)
or
"compact" spline
The behaviour of the bicubic Hermite interpolant is analogous to the 1d
case. The required input information is the function value and
derivatives at each grid point:
x(i), i = 1 to Nx, x(i).lt.x(i+1) -- the x grid
y(j), j = 1 to Ny, y(j).lt.y(j+1) -- the y grid
f(0,i,j) = f(x(i),y(j))
f(1,i,j) = [df/dx](x(i),y(j))
f(2,i,j) = [df/dy](x(i),y(j))
f(3,i,j) = [d2f/dxdy](x(i),y(j))
Bicubic Hermite evaluation at (x,y) is done by direct construction,
as follows. Suppose that (x,y) is in the grid cell whose lower left
corner is (x(i),y(j)).
Let hx(i)=(x(i+1)-x(i)),
hxi(i)=1/hx(i),
hy(j)=(y(j+1)-y(j)),
hyi(j)=1/hy(j).
x~=(x-x(i))*hxi(i)
y~=(y-y(j))*hyi(j)
The normalized parameters (x~,y~) describe the position of (x,y) within
the grid cell.
The Hermite basis functions alpha, alphabar, beta, and betabar
(defined in the previous section) are then used:
Let
ax=alpha(x~)
axbar=alphabar(x~)
bx=beta(x~)
bxbar=betabar(x~)
ay=alpha(y~)
aybar=alphabar(y~)
by=beta(y~)
bybar=betabar(y~)
then
h(x,y) = aybar*(axbar*f(0,i,j)+ax*f(0,i+1,j))
+ ay*(axbar*f(0,i,j+1)+ax*f(0,i+1,j+1))
+ hx(i)*(
aybar*(bxbar*f(1,i,j)+bx*f(1,i+1,j))
+ ay*(bxbar*f(1,i,j+1)+bx*f(1,i+1,j+1)) )
+ hy(j)*(
bybar*(axbar*f(2,i,j)+ax*f(2,i+1,j))
+ by*(axbar*f(2,i,j+1)+ax*f(2,i+1,j+1)) )
+ hx(i)*hy(j)*
bybar*(bxbar*f(3,i,j)+bx*f(3,i+1,j))
+ by*(bxbar*f(3,i,j+1)+bx*f(3,i+1,j+1)) )
with similar expressions, involving derivatives axbar', ax',
bxbar', and bx', for d[h(x,y)]/dx, and involving derivatives aybar',
ay', bybar', and by', for d[h(x,y)]/dy, and all polynomial
derivatives for d2[h(x,y)]/dxdy. Full details can be seen
in the PSPLINE library source code: herm2ev.for.
From the definitions of the basis functions and their it follows
that all the necessary continuity and differentiability conditions
are met. For example, h(x(i),y(j))=f(0,i,j) (as required) because
x(i) -> x~=0, y(i) -> y~=0, => axbar=aybar=1, ax=ay=0, and
bx=bxbar=by=bybar=0.
Similarly, d[h(x(i),y(j+1)]/dy = f(2,i,j+1) (as required) because
x~=0, y~=1 and the known simple behaviour of the basis functions and
their derivatives.
Inside each cell, h(x,y) is locally a bicubic polynomial patch, and
is multiply continuously differentiable.
Along any cell boundary, the formulae for h(x,y) and its 1st derivatives
reduce to expressions which depend only on the data at the endpoints
of the line segment defining the boundary. These limiting expressions
are identical for cells on either side of the boundary segment. This is
what yields the continuity of h(x,y) and all its 1st derivatives across
these boundaries and hence over the entire gridded region.
The PSPLINE library provides routines for evaluation of 2d Hermite
interpolants and their first derivatives -- see herm2ev.for.
Hermite interpolation is different than spline interpolation in
that (a) no global matrix solution is required for evaluation of
coefficients, and (b) continuity properties are maintained for ANY
choice of derivatives at the grid points.
Thus, if derivative information is not available analytically, it
may be provided by other suitable (i.e. numerical) means.
The PSPLINE library Hermite software includes a choice of routines
for numerical evaluation of derivatives-- see the section on software.
The construction method using 1d Hermite cubic basis functions, described
for 2d bicubic Hermite, generalizes in a straightforward manner to the 3d
case as well. The details will not be spelled out here -- they can be
seen in the PSPLINE library source code: herm3ev.for.
The input data needed for tricubic Hermite evaluation are:
x(i), i = 1 to Nx, x(i).lt.x(i+1) -- the x grid
y(j), j = 1 to Ny, y(j).lt.y(j+1) -- the y grid
z(k), k = 1 to Nz, z(k).lt.z(k+1) -- the y grid
f(0,i,j,k) = f(x(i),y(j),z(k))
f(1,i,j,k) = [df/dx](x(i),y(j),z(k))
f(2,i,j,k) = [df/dy](x(i),y(j),z(k))
f(3,i,j,k) = [df/dz](x(i),y(j),z(k))
f(4,i,j,k) = [d2f/dxdy](x(i),y(j),z(k))
f(5,i,j,k) = [d2f/dxdz](x(i),y(j),z(k))
f(6,i,j,k) = [d2f/dydz](x(i),y(j),z(k))
f(7,i,j,k) = [d3f/dxdydz](x(i),y(j),z(k))
The PSPLINE library provides routines for evaluation of 3d Hermite
interpolants and their first derivatives -- see herm3ev.for.
Hermite interpolation is different than spline interpolation in
that (a) no global matrix solution is required for evaluation of
coefficients, and (b) continuity properties are maintained for ANY
choice of derivatives at the grid points.
Thus, if derivative information is not available analytically, it
may be provided by other suitable (i.e. numerical) means.
The PSPLINE library Hermite software includes a choice of routines
for numerical evaluation of derivatives-- see the section on software.
Interpolation software often represents a performance critical
component of real applications. Here, in decreasing order of
importance, are some tips for obtaining the best possible
performance from pspline:
* where possible, do "vectorized" spline evaluations
(including zone lookup), i.e. use routines which carry out
multiple spline evaluations in a single call.
"Zone lookup" plays a critical role in the performance of
interpolation systems.
* where possible, use evenly spaced grids. Otherwise...
* using custom zone lookup software may help when special
properties of the grid are known in advance. Otherwise...
* use "(r8)genxpkg" to create a grid data object for each
grid dimension, and use the corresponding "genxpkg-
enabled" set of evaluation routines.
-> On all workstations tested so far, genxpkg algorithm #3,
piecewise-linear indexing function, is the fastest, out-
performing traditional binary search methods by almost a
factor of two.
-> if the target point locations for interpolation are
slowly varying, selecting the genxpkg option to "use
previous result" can be a win (~30% performance gain
on typical workstations). On the other hand, if the
target point locations vary rapidly, this option will
exact a small penalty (~5% performance loss on typical
workstations). Caution: in an evaluation loop, "use
previous result" creates a dependency on the previous
loop iteration, which may break vectorization and
exact a heavy performance penalty, on some machines.
There is a trade-off between memory and speed in this software...
* use "explicit" splines rather than "compact" splines, when
the memory cost is not too high. Prefer "explicit" splines
to Hermite methods, unless the spline setup cost is a
significant factor, or data noise is a consideration.
The following step will involve too much labor, except where the
interpolation performance is extraordinarily critical
* for the ultimate in performance, hand inline the spline
evaluation computation into your application (but note this
gets tricky, especially at higher dimensionalities)!
When moving to a finer grid, C2 cubic splines will converge much
more rapidly than C1 Hermite to an accurate reproduction of the
underlying model function, if the function is smooth and the
representative data points chosen sample the model function
with machine precision accuracy.
But: when fitting "noisy experimental data", splines can cause
amplification and spreading of the noise, called "ringing". The
same comment applies if the underlying model function is not
really C2 on the resolution of the given numerical grid. If
noisy data is expected, Hermite methods should be given preference.
The Akima Hermite method, in particular, has well demonstrated
variation minimizing properties.
[Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1]
Most users should follow this strategy:
* use EZspline
But, if the f77 interface must be used, then,
* use genxpkg to control zone lookup algorithms and precompute
grid-specific data which is useful to the lookup/evaluation
software, and, use vectorized spline lookup/evaluation routines,
or...
* provide your own vectorized zone lookup algorithm, and use the
vectorized spline evaluation routines.
fortran-90 interface with good performance and ease of use features
by Alex Pletzer...
for details, see the separate EZspline help documents. Links to
these documents can be found at the pspline web page:
http://w3.pppl.gov/NTCC/PSPLINE/pspline.html
The subtopics describe the pspline "standard routines" for Spline
and Hermite interpolation. See also the "EZspline_Software"
section for the description of Alex Pletzer's f90 interface to
pspline, which will probably be easier to use for new code.
Note: the routines described are coded for single precision
(floating point arguments and arrays are declared REAL ...),
but...
The PSPLINE library also contains a complete set of double
precision routines (floating point arguments and arrays are
declared REAL*8 ...). These routines have the same names as
the single precision routines except that the prefix "R8"
prepended.
For example:
standard name REAL*8 name
one-d spline setup routine CSPLINE R8CSPLINE
two-d spline setup routine BCSPLINE R8BCSPLINE
The declaration "REAL*8" instead of "DOUBLE PRECISION" is used
for reasons of portability to typical Cray environments.
The pspline software includes a "public" module with a definition
of all the (non-EZspline) pspline subroutine calls. Fortran-90
users can add the line
use pspline_calls
to their source, for compile-time checking of argument lists
of pspline subroutine calls.
Caution: as with any f90 interface module, some calls that are
functionally correct may be flagged. This happens e.g. if an
element of an array is passed to a dummy argument as an "address"
indicating that that array element and subsequent elements in the
caller's array contain the necessary data. Unfortunately, f90
sees this as an attempt to pass a scalar where an array is expected.
Note that the use of the pspline_calls modules is entirely optional:
it is essentially a debugging aid.
This routine takes a grid array x(1:nx) and prepares a dataset
xpkg(1:nx,4) imbedding precomputed data specific to the x
axis, and option specifications for pspline zone lookup software.
Usage:
real x(nx) ! the grid (in) -- strict ascending
real xpkg(nx,4) ! the "package" (out)
integer iper ! periodicity flag (in)
integer imsg ! flag for out-of-range warnings (in)
integer itol ! range tolerance option (in)
real ztol ! range tolerance, if itol is set (in)
integer ialg ! algorithm selection code (in)
integer ier ! completion code, 0=OK (out)
call genxpkg(nx,x,xpkg,iper,imsg,itol,ztol,ialg,ier)
c!--------------------
c! set iper=1 if x is a periodic coordinate, e.g. poloidal or
c! toroidal angle coordinate in a toroidal system.
c! if x is periodic, a mod function is applied to
c! normalize input x locations that are out of the
c! given range. x(1:nx) is presumed to cover 1
c! whole period.
c! set iper=0 if x is not a periodic coordinate
c!
c! recommendation: set correctly for the given x coordinate
c!--------------------
c! set imsg=1 for error messages when a given input x location
c! is beyond the range of x(1:nx).
c! set imsg=0 for no error message, but, regardless of imsg, a
c! warning flag will be set by the lookup routine
c!
c! recommendation: set imsg=0 but check warning flags in
c! routines which call for zone lookup / spline
c! evaluation.
c!--------------------
c! set itol=1 to specify relative tolerance for determining
c! even spacing, and for out-of-range errors
c! set itol=0 to use the default relative tolerance, (5.0e-7).
c!
c! set ztol= <the relative tolerance> if itol=1
c!
c! absolute tolerance is then
c! ztola = ztol*max(abs(x(1)),abs(x(nx)));
c! a data point within the range x(1)-ztola to x(nx)+ztola
c! is considered to be "in range".
c!
c! recommendation: the default is OK for most purposes.
c!--------------------
c! to set the search algorithm, in case the grid is unevenly
c! spaced:
c!
c! set ialg= +/- 1: <1/h(j)> "Newton" method
c! +/- 2: binary search
c! +/- 3: piecewise linear indexing function
c!
c! choose the negative value to allow use of the results of
c! the search for the preceding target point, as a start point
c! for the search for the subsequent target point.
c!
c! recommendation: for workstations: ialg = -3 ... is fastest.
c! if the grid is determined to be evenly spaced, this fact
c! will be recorded in xpkg, and the lookup software will
c! take advantage of this, regardless of the setting of "ialg".
c!--------------------
pspline comes with test programs "lookup_test" and "r8lookup_test"
to exercise the even-spaced and uneven-spaced grid zone lookup
algorithms, measuring cpu time consumption.
Tests were run (12 Apr 2000) on a COMPAQ/DEC Alpha 500au using
the DEC FORTRAN optimizing compiler (version 5.1) under
DIGITAL UNIX v4.0.
(a) "vectorization" has a strong performance impact: for example,
(500000 lookups per test)
...even spacing: NO
...algorithm: linear indexing fcn
...DO NOT use previous search result.
%genpkg(ngrid,grid,gpkg, 1,1,0,zdum, 3,ier): ier=
0
$$$$$> cpu time (secs) = 1.312721 vector size = 1
$$$$$> cpu time (secs) = 0.2918243 vector size = 10
$$$$$> cpu time (secs) = 0.1756783 vector size = 100
$$$$$> cpu time (secs) = 0.1668968 vector size = 1000
(this compares 500000 lookups with varying levels of vectorization,
ranging from 500000 single lookup calls to 500 calls with 1000
lookups each).
(b) for vector size = 10, testing 500000 lookups with on an
unevenly spaced grid. Here, numbers in parentheses are from a
Linux P-II 300MHz system running Red Hat 6.1 and using the fujitsu
fortran compiler v1.0.
cpu time, seconds
slowly varying target vector:
-----------------------------
algorithm init search from each iteration
previous iteration independent
Pseudo-Newton .304 (.73) .469 (1.07)
Binary-Search .381 (.85) .508 (.99)
Index-Function .208 (.49) .297 (.67)
rapidly varying target vector:
------------------------------
algorithm init search from each iteration
previous iteration independent
Pseudo-Newton .534 (1.19) .515 (1.20)
Binary-Search .566 (1.14) .537 (1.06)
Index-Function .321 (.71) .286 (.67)
(date of test: 12-Apr 2000)
Pseudo-Newton (genxpkg ialg = +/- 1) -- start in zone j. If that
is not the right zone, use delta(x)/h(j) to estimate which
zone to look in next. delta(x) = distance of target point
from zone j boundary; h(j) = x(j+1)-x(j).
Binary-Search (genxpkg ialg = +/- 2) -- start in zone nx/2. If
that is not the right zone, then, if the zone is too low try
3*nx/4 else try 1*nx/4. Keep subdividing the possible set of
zones by a factor of 2 per step, until the right zone is reached.
We expect log[2](nx) performance scaling of this algorithm.
Index-Function (genxpkg ialg = +/- 3) -- we have pretabulated
on an evenly spaced grid covering x(1) to x(nx) the index
(integer part) and offset with/in zone (fractional part) into
the original, unevenly spaced grid. To do a lookup, we carry out
a piecewise linear interpolation of this indexing function and
extract the integer part of the result. If this yields the
correct zone in the original grid, we are done. If not, we
correct by shifting right or left one zone at a time.
This routine takes a vector of target locations and determines
the corresponding vectors of zone indices and displacements.
It uses "xpkg" (output of subroutine genxpkg-- see above) to
acquire the grid data and lookup algorithm options.
The "imode" switch controls whether to generate location
information required by the explicit spline evaluation
routines (imode=1) or in the form required by compact
spline, Hermite, and piecewise linear routines (imode=2).
Usage:
**call genxpkg first**
real xvec(ivec) ! target vector (in)
real xpkg(nx,4) ! grid xpkg (in)
integer imode ! mode switch (in)
integer iv(ivec) ! list of zone indices (out)
real dxn(ivec) ! displacements w/in zones (out)
real hv(ivec) ! zone widths (out)
real hiv(ivec) ! inverse zone widths (out)
integer iwarn ! warning flag, 0=normal (out)
call xlookup(ivec,xvec,nx,xpkg,imode,iv,dxn,hv,hiv,iwarn)
c! xpkg: must be set up by prior call to genxpkg.
c! xpkg(1:nx,1) = x(1:nx), where x(1:nx) is the original
c! grid.
c!
c! imode switch:
c!
c! if imode=1, then, on output:
c! iv(i) contains the zone in the original grid x(1:nx)
c! in which xvec(i) is contained;
c! i.e. x(i).le.xvec(i).le.x(i+1)
c! iv(i) is always in the range [1,nx-1].
c!
c! dxn(i) = xvec(i) - x(i)
c!
c! (hv(i) and hiv(i) are never set)
c!
c! *** use this for non-compact spline evaluation ONLY ***
c!
c! if imode=2, then, on output:
c! iv(i) contains the zone in the original grid x(1:nx)
c! in which xvec(i) is contained;
c! i.e. x(i).le.xvec(i).le.x(i+1)
c! iv(i) is always in the range [1,nx-1].
c!
c! dxn(i) = (xvec(i) - x(i))/hv(i)
c!
c! hv(i) = (x(i+1)-x(i))
c!
c! hiv(i) = 1/hv(i)
c!
c! *** use this for everything except non-compact splines ***
c!
c! out-of-range warnings (iwarn):
c!
c! iwarn=0 on output: no warnings
c! iwarn=M on output: M points in xvec(1:ivec) were OUT OF RANGE
c!
c! If x(1:nx) is specified periodic (cf genxpkg call), range errors
c! cannot occur; a mod function is used instead to normalize the
c! range of the periodic coordinate; iwarn=0 on output.
c!
c! points xvec(j) out of range low are treated as if xvec(j)=x(1).
c! points xvec(k) out of range high are treated as if xvec(k)=x(nx).
c!
c! If out of range points are within a tolerance of the boundary,
c! the warning flag is not set. The tolerance was set by the
c! genxpkg call.
c!
Routines for setting up spline coefficients:
compact splines memory requirement
---------------------------------------------
1d: mkspline 2*Nx
2d: mkbicub 4*Nx*Ny
3d: mktricub 8*Nx*Ny*Nz
==> these routines create the spline coefficients and any other
information needed by the lookup / evaluation routines --
with various boundary condition control options.
Vectorized Evaluation Routines (recommended):
=============================================
given a location vector, do lookup & evaluation:
vecspline vecbicub vectricub
given a grid (different than the original grid), to
map the function from the old grid to the new grid:
gridspline gridbicub gridtricub
if zone location data is already computed, to do just
the vectorized spline evaluation:
fvspline fvbicub fvtricub
Scalar Routines:
================
Routines for compact spline function evaluation (scalar, grid
lookup included):
evspline evbicub evtricub
subroutine mkspline(x,nx,fspl,ibcxmin,bcxmin,ibcxmax,bcxmax,
> ilinx,ier)
C
C make a 2-coefficient 1d spline
C
C only 2 coefficients, the data and its 2nd derivative, are needed to
C fully specify a spline. See e.g. Numerical Recipies in Fortran-77
C (2nd edition) chapter 3, section on cubic splines.
C
C input:
integer nx ! no. of data points
real x(nx) ! x axis data, strict ascending order
C
C input/output:
real fspl(2,nx) ! f(1,*): data in; f(2,*): coeffs out
C
C f(1,j) = f(x(j)) on input (unchanged on output)
C f(2,j) = f''(x(j)) (of interpolating spline) (on output).
C
C ...boundary conditions...
C
C input:
C
integer ibcxmin ! b.c. flag @ x=xmin=x(1)
real bcxmin ! b.c. data @xmin
C
integer ibcxmax ! b.c. flag @ x=xmax=x(nx)
real bcxcmax ! b.c. data @xmax
C
C ibcxmin=-1 -- periodic boundary condition
C (bcxmin,ibcxmax,bcxmax are ignored)
C
C the output spline s satisfies
C s'(x(1))=s'(x(nx)) ..and.. s''(x(1))=s''(x(nx))
C
C if non-periodic boundary conditions are used, then the xmin and xmax
C boundary conditions can be specified independently:
C
C ibcxmin (ibcxmax) = 0 -- this specifies a "not a knot" boundary
C condition, see "cubsplb.for". This is a common way
C for inferring a "good" spline boundary condition
C automatically from data in the vicinity of the
C boundary. ... bcxmin (bcxmax) are ignored.
C
C ibcxmin (ibcxmax) = 1 -- boundary condition is to have s'(x(1))
C ( s'(x(nx)) ) match the passed value bcxmin (bcxmax).
C
C ibcxmin (ibcxmax) = 2 -- boundary condition is to have s''(x(1))
C ( s''(x(nx)) ) match the passed value bcxmin (bcxmax).
C
C ibcxmin (ibcxmax) = 3 -- boundary condition is to have s'(x(1))=0.0
C ( s'(x(nx))=0.0 )
C
C ibcxmin (ibcxmax) = 4 -- boundary condition is to have s''(x(1))=0.0
C ( s''(x(nx))=0.0 )
C
C ibcxmin (ibcxmax) = 5 -- boundary condition is to have s'(x(1))
C ( s'(x(nx)) ) match the 1st divided difference
C e.g. at x(1): d(1)/h(1), where
C d(j)=f(1,j+1)-f(1,j)
C h(j)=x(j+1)-x(j)
C
C ibcxmin (ibcxmax) = 6 -- BC is to have s''(x(1)) ( s''(x(nx)) )
C match the 2nd divided difference
C e.g. at x(1):
C e(1) = [d(2)/h(2) - d(1)/h(1)]/(0.5*(h(1)+h(2)))
C
C ibcxmin (ibcxmax) = 7 -- BC is to have s'''(x(1)) ( s'''(x(nx)) )
C match the 3rd divided difference
C e.g. at x(1): [e(2)-e(1)]/(0.33333*(h(1)+h(2)+h(3)))
C
C output:
C
integer ilinx ! =1: hint, x axis is ~evenly spaced
C
C let dx[avg] = (x(nx)-x(1))/(nx-1)
C let dx[j] = x(j+1)-x(j), for all j satisfying 1.le.j.lt.nx
C
C if for all such j, abs(dx[j]-dx[avg]).le.(1.0e-3*dx[avg]) then
C ilinx=1 is returned, indicating the data is (at least nearly)
C evenly spaced. Even spacing is useful, for speed of zone lookup,
C when evaluating a spline.
C
C if the even spacing condition is not satisfied, ilinx=2 is returned.
C
integer ier ! exit code, 0=OK
C
C an error code is returned if the x axis is not strict ascending,
C or if nx.lt.4, or if an invalid boundary condition specification was
C input.
subroutine mkbicub(x,nx,y,ny,f,nf2,
> ibcxmin,bcxmin,ibcxmax,bcxmax,
> ibcymin,bcymin,ibcymax,bcymax,
> ilinx,ilinth,ier)
From mkbicub.for:
subroutine mkbicub(x,nx,y,ny,f,nf2,
> ibcxmin,bcxmin,ibcxmax,bcxmax,
> ibcymin,bcymin,ibcymax,bcymax,
> ilinx,iliny,ier)
C
C setup a bicubic spline; store coefficients in compact form
C (as per suggestion of L. Zakharov, PPPL, Feb. 1999)
C 4 coeffs per grid point: f,fxx,fyy,fxxyy
C
C
C input:
integer nx ! length of x vector
integer ny ! length of y vector
real x(nx) ! x vector, strict ascending
real y(ny) ! y vector, strict ascending
C
integer nf2 ! 2nd dimension of f, nf2.ge.nx
C input/output:
C
real f(4,nf2,ny) ! data & spline coefficients
C
C on input: f(1,i,j) = f(x(i),y(j))
C on output: f(1,i,j) unchanged
C f(2,i,j) = d2f/dx2(x(i),y(j))
C f(3,i,j) = d2f/dy2(x(i),y(j))
C f(4,i,j) = d4f/dx2dy2(x(i),y(j))
C
C and the interpolation formula for (x,y) in
C (x(i),x(i+1))^(y(j),y(j+1)) is:
C hx = x(i+1)-x(i) hy = y(j+1)-y(j)
C dxp= (x-x(i))/hx dxm= 1-dxp dxp,dxm in (0,1)
C dyp= (x-x(i))/hx dym= 1-dyp dyp,dym in (0,1)
C dx3p = dxp**3-dxp dx3m = dxm**3-dxm dxp3,dxm3 in (0,1)
C
C finterp = dxm*(dym*f(1,i,j)+dyp*f(1,i,j+1))
C +dxp*(dym*f(1,i+1,j)+dyp*f(1,i+1,j+1))
C +1/6*hx**2*
C dx3m*(dym*f(2,i,j)+dyp*f(2,i,j+1))
C +dx3p*(dym*f(2,i+1,j)+dyp*f(2,i+1,j+1))
C +1/6*hy**2*
C dxm*(dy3m*f(3,i,j)+dy3p*f(3,i,j+1))
C +dxp*(dy3m*f(3,i+1,j)+dy3p*f(3,i+1,j+1))
C +1/36*hx**2*hy**2*
C dx3m*(dym*f(4,i,j)+dyp*f(4,i,j+1))
C +dx3p*(dym*f(4,i+1,j)+dyp*f(4,i+1,j+1))
C
C where the f(2:4,*,*) are cleverly evaluated to assure
C (a) finterp is continuous and twice differentiable across all
C grid cell boundaries, and
C (b) all boundary conditions are satisfied.
C
C the boundary conditions specification options are:
C inputs:
C
integer ibcxmin ! bc flag for x=xmin
real bcxmin(ny) ! bc data vs. y at x=xmin
integer ibcxmax ! bc flag for x=xmax
real bcxmax(ny) ! bc data vs. y at x=xmax
C
integer ibcymin ! bc flag for y=ymin
real bcymin(nx) ! bc data vs. x at y=ymin
integer ibcymax ! bc flag for y=ymax
real bcymax(nx) ! bc data vs. x at y=ymax
C
C with interpretation:
c ibcxmin -- indicator for boundary condition at x(1):
c bcxmin(...) -- boundary condition data
c =-1 -- periodic boundary condition
c =0 -- use "not a knot"
c =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
c =2 -- match 2nd deriv., specified at x(1),th(ith) by bcxmin(ith)
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
c =5 -- match 1st derivative to 1st divided difference
c =6 -- match 2nd derivative to 2nd divided difference
c =7 -- match 3rd derivative to 3rd divided difference
c (for more detailed definition of BCs 5-7, see the
c comments of subroutine mkspline)
c NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
c
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c NOTE: if ibcxmin=-1, ibcxmax is ignored! ...the BC is periodic.
c
C and analogous interpretation for ibcymin,bcymin,ibcymax,bcymax
C (df/dy or d2f/dy2 boundary conditions at y=ymin and y=ymax).
C
C
C and output arguments
C
integer ilinx ! =1: x grid is "nearly" equally spaced
integer iliny ! =1: y grid is "nearly" equally spaced
C
C ilinx and iliny are set to zero if corresponding grids are not equally
C spaced
C
C and completion code
C
integer ier ! =0 on exit if there is no error.
C
C if there is an error, ier is set and a message is output on unit 6.
C these are considered programming errors in the calling routine.
C
C possible errors:
C x(...) not strict ascending
C y(...) not strict ascending
C nx .lt. 4
C ny .lt. 4
C invalid boundary condition flag
C workspace too small
subroutine mktricub(x,nx,y,ny,z,nz,f,nf2,nf3,
> ibcxmin,bcxmin,ibcxmax,bcxmax,inb1x,
> ibcymin,bcymin,ibcymax,bcymax,inb1y,
> ibczmin,bczmin,ibczmax,bczmax,inb1z,
> ilinx,iliny,ilinz,ier)
c
c setup a tricubic spline; store coefficients in compatct form
c (as per suggestion of L. Zakharov, PPPL, Feb. 1999)
C 8 coeffs per (x,y,z) grid point:
C f,fxx,fyy,fzz,fxxyy,fxxzz,fyyzz,fxxyyzz
C
C input:
integer nx ! length of x vector
integer ny ! length of y vector
integer nz ! length of z vector
real x(nx) ! x vector, strict ascending
real y(ny) ! y vector, strict ascending
real z(nz) ! z vector, strict ascending
c
integer nf2 ! 2nd dim. of f array, nf2.ge.nx
integer nf3 ! 3rd dim. of f array, nf3.ge.ny
c
c input/output:
c
real f(8,nf2,nf3,nz) ! data and spline coefficients
c
C on input: f(1,i,j,k) = f(x(i),y(j),z(k))
C on output: f(1,i,j,k) unchanged
C f(2,i,j,k) = d2f/dx2(x(i),y(j),z(k))
C f(3,i,j,k) = d2f/dy2(x(i),y(j),z(k))
C f(4,i,j,k) = d2f/dz2(x(i),y(j),z(k))
C f(5,i,j,k) = d4f/dx2dy2(x(i),y(j),z(k))
C f(6,i,j,k) = d4f/dx2dz2(x(i),y(j),z(k))
C f(7,i,j,k) = d4f/dy2dz2(x(i),y(j),z(k))
C f(8,i,j,k) = d6f/dx2dy2dz2(x(i),y(j),z(k))
C
C there is a rather Hermite like interpolation formula to go with
C this-- see evtricub.for. Also the bicubic formula is given in
C mkbicubw.for; the tricubic formula is precisely analogous.
C
C boundary condition data
C inputs:
integer inb1x ! 1st dim of xmin & xmax bc arrays
integer inb1y ! 1st dim of ymin & ymax bc arrays
integer inb1z ! 1st dim of zmin & zmax bc arrays
C
integer ibcxmin,ibcxmax ! BC type flag @xmin, xmax
integer ibcymin,ibcymax ! BC type flag @ymin, ymax
integer ibczmin,ibczmax ! BC type flag @zmin, zmax
C
real bcxmin(inb1x,nz)
real bcxmax(inb1x,nz) ! xmin & xmax BC data, ny x nz
real bcymin(inb1y,nz)
real bcymax(inb1y,nz) ! ymin & ymax BC data, nx x nz
real bczmin(inb1z,ny)
real bczmax(inb1z,ny) ! zmin & zmax BC data, nx x ny
c
c where BC data is not required, dummy scalars may be passed.
C the ibc* flags determine whether BC data isneeded.
c
c BC data: bcxmin & bcxmax: BC vs. y,z @xmin,xmax
C bcymin & bcymax: BC vs. x,z @ymin,ymax
C bczmin & bczmax: BC vs. x,y @zmin,zmax
C
c ibcxmin -- indicator for boundary condition at xmin=x(1):
c bcxmin(...) -- boundary condition data
c =-1 -- use periodic boundary condition
c =0 -- use "not a knot"
c =1 -- match slope, specified at x(1),y(iy),z(iz) by bcxmin(iy,iz)
c =2 -- match 2nd derivative, specified at x(1),y(iy),z(iz)
c by bcxmin(iy,iz
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all y(j)
c =4 -- boundary condition is d2f/dx2=0 at x(1), all y(j)
c =5 -- match 1st derivative to 1st divided difference
c =6 -- match 2nd derivative to 2nd divided difference
c =7 -- match 3rd derivative to 3rd divided difference
c (for more detailed definition of BCs 5-7, see the
c comments of subroutine mkspline)
c ***NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
c
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c NOTE: if ibcxmin=-1 then the periodic BC applies on both sides
c and ibcxmax, bcxmax are ignored.
c inb1x -- 1st dimension of bcxmin, bcxmax: if ibcxmin or ibcxmax
c .gt. 0 this must be .ge. ny.
c
c interpretation of ibcymin,bcymin,ibcymax,bcymax,inb1y
c is same as with ibcxmin,... ** inb1y.ge.nx required **
c
c interpretation of ibczmin,bczmin,ibczmax,bczmax,inb1z
c is same as with ibcxmin,... ** inb1z.ge.nx required **
c
c the explicit bdy condition arrays are referenced only if the
c corresponding "ibc" flag values are set to 1 or 2.
c
c output:
integer ilinx ! x vector equal spacing flag
integer iliny ! y vector equal spacing flag
integer ilinz ! z vector equal spacing flag
c
c ilinx -- =1 on output if x(nx) nearly evenly spaced (tol=1e-3)
c iliny -- =1 on output if y(ny) evenly spaced (tol=1e-3)
c ilinz -- =1 on output if z(nz) evenly spaced (tol=1e-3)
c
integer ier ! exit code
c ier -- completion code, 0 for normal
subroutine vecspline(ict,ivec,xvec,ivd,fval,nx,xpkg,fspl,
> iwarn,ier)
c
c vectorized spline evaluation routine -- 1d *compact* spline
c 1. call vectorized zone lookup routine
c 2. call vectorized spline evaluation routine
c
c--------------------------
c input:
integer ict(3) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for d2f/dx2 (don't evaluate if ict(3)=0)
c
c set ict(1)=3 to get d3f/dx3 (only)
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c
real xvec(ivec) ! x-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx ! dimension of spline x grid
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real fspl(2,nx) ! (compact) spline coefficients
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
subroutine vecbicub(ict,ivec,xvec,yvec,ivd,fval,
> nx,xpkg,ny,ypkg,fspl,inf2,
> iwarn,ier)
c
c vectorized spline evaluation routine -- 2d *compact* spline
c 1. call vectorized zone lookup routine
c 2. call vectorized spline evaluation routine
c
c--------------------------
c input:
integer ict(6) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
c ict(4)=1 for d2f/dx2 (don't evaluate if ict(4)=0)
c ict(5)=1 for d2f/dy2 (don't evaluate if ict(5)=0)
c ict(6)=1 for d2f/dxdy (don't evaluate if ict(6)=0)
c the number of non zero values ict(1:6)
c determines the number of vector outputs...
c
c new dmc December 2005 -- access to higher derivatives (even if not
c continuous-- but can only go up to 3rd derivatives on any one coordinate.
c if ict(1)=3 -- want 3rd derivatives
c ict(2)=1 for d3f/dx3
c ict(3)=1 for d3f/dx2dy
c ict(4)=1 for d3f/dxdy2
c ict(5)=1 for d3f/dy3
c number of non-zero values ict(2:5) gives no. of outputs
c if ict(1)=4 -- want 4th derivatives
c ict(2)=1 for d4f/dx3dy
c ict(3)=1 for d4f/dx2dy2
c ict(4)=1 for d4f/dxdy3
c number of non-zero values ict(2:4) gives no. of outputs
c if ict(1)=5 -- want 5th derivatives
c ict(2)=1 for d5f/dx3dy2
c ict(3)=1 for d5f/dx2dy3
c number of non-zero values ict(2:3) gives no. of outputs
c if ict(1)=6 -- want 6th derivatives
c d6f/dx3dy3 -- one value is returned.
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c
c list of (x,y) pairs:
c
real xvec(ivec) ! x-locations at which to evaluate
real yvec(ivec) ! y-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx,ny ! dimension of spline grids
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real ypkg(ny,4) ! y grid "package" (cf genxpkg)
integer inf2 ! fspl 3rd array dimension, .ge.nx
real fspl(0:3,inf2,ny) ! (compact) spline coefficients
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
c
subroutine vectricub(ict,ivec,xvec,yvec,zvec,ivd,fval,
> nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf2,inf3,
> iwarn,ier)
c
c vectorized spline evaluation routine -- 3d *compact* spline
c 1. call vectorized zone lookup routine
c 2. call vectorized spline evaluation routine
c
c--------------------------
c input:
integer ict(10) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for df/dy (don't evaluate if ict(3)=0)
c ict(4)=1 for df/dy (don't evaluate if ict(4)=0)
c ict(5)=1 for d2f/dx2 (don't evaluate if ict(5)=0)
c ict(6)=1 for d2f/dy2 (don't evaluate if ict(6)=0)
c ict(7)=1 for d2f/dz2 (don't evaluate if ict(7)=0)
c ict(8)=1 for d2f/dxdy (don't evaluate if ict(8)=0)
c ict(9)=1 for d2f/dxdz (don't evaluate if ict(9)=0)
c ict(10)=1 for d2f/dydz (don't evaluate if ict(10)=0)
c
c (new dmc Dec 2005 -- higher derivatives available)
c ict(1)=3 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:8) select: fxxy fxxz fxyy fxyz fxzz fyyz fyzz
c ->note ict(1)=3, ict(5)=1 gives fxyz = d3f/dxdydz
c ict(1)=-3 --> 3rd derivative, 3 in one coordinate
c ict(2:4) select: fxxx fyyy fzzz
c ict(1)=4 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:7) select: fxxyy fxxyz fxxzz fxyyz fxyzz fyyzz
c ict(1)=-4 --> 3rd derivative, 3 in one coordinate
c ict(2:7) select: fxxxy fxxxz fxyyy fxzzz fyyyz fyzzz
c ict(1)=5 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:4) select: fxxyyz fxxyzz fxyyzz
c ict(1)=-5 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyy fxxxyz fxxxzz fxxyyy fxxzzz
c fxyyyz fxyzzz fyyyzz fzzzyy
c ict(1)=6 --> 3rd derivative, .le.2 diff. in any coordinate
c fxxyyzz
c ict(1)=-6 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyyy fxxxyyz fxxxyzz fxxxyyz
c fxxyyyz fxxyzzz fxyyyzz fxyyzzz fyyyzzz
c ict(1)=-7 --> 7th derivative
c ict(2:7) select: fxxxyyyz fxxxyyzz fxxxyzzz
c fxxyyyzz fxxyyzzz fxyyyzzz
c ict(1)=-8 --> 8th derivative
c ict(2:4) select: fxxxyyyzz fxxxyyzzz fxxyyyzzz
c ict(1)=-9 --> 9th derivative: fxxxyyyzzz
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c
c list of (x,y,z) triples:
c
real xvec(ivec) ! x-locations at which to evaluate
real yvec(ivec) ! y-locations at which to evaluate
real zvec(ivec) ! z-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values as per 1st non-zero ict(...) element
c fval(1:ivec,2) -- values as per 2nd non-zero ict(...) element
c --etc--
c
c input:
integer nx,ny,nz ! dimension of spline grids
real xpkg(nx,4) ! x grid "package" (cf genxpkg)
real ypkg(ny,4) ! y grid "package" (cf genxpkg)
real zpkg(nz,4) ! z grid "package" (cf genxpkg)
integer inf2 ! fspl 4th array dimension, .ge.nx
integer inf3 ! fspl 5th array dimension, .ge.ny
real fspl(0:7,inf2,inf3,nz) ! (compact) spline coefficients
c
c output:
c condition codes, 0 = normal return
integer iwarn ! =1 if an x value was out of range
integer ier ! =1 if argument error detected
subroutine gridspline(x_newgrid,nx_new,f_new,nx,xpkg,fspl,
> iwarn,ier)
c
c regrid a spline function f defined vs. x as in xpkg
c to a new grid, given by x_newgrid.
c
c set warning flag if the range x_newgrid exceeds the range of the
c original xpkg.
c
c (xpkg -- see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new grid
c
c output:
c
real f_new(nx_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
real fspl(nx,2) ! compact spline coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
c
c--------------------------------------------
c local
c
integer ict(3)
c
data ict/1,0,0/
c
c--------------------------------------------
c
call vecspline(ict,nx_new,x_newgrid,nx_new,f_new,nx,xpkg,fspl,
> iwarn,ier)
c
return
end
subroutine gridbicub(
> x_newgrid,nx_new,y_newgrid,ny_new,f_new,if1,
> nx,xpkg,ny,ypkg,fspl,inf3,iwarn,ier)
c
c regrid a spline function f defined vs. x,y as in xpkg,ypkg
c to a new grid, given by x_newgrid, y_newgrid
c
c set warning flag if the range x_newgrid, y_newgrid
c exceeds the range of the original xpkg,ypkg.
c
c (xpkg,ypkg -- axis data, see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new x grid
real y_newgrid(ny_new) ! new y grid
c
c output:
c
integer if1 ! 1st dimension of f_new
real f_new(if1,ny_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
integer ny ! size of old grid
real ypkg(ny,4) ! old grid "package"
c
integer inf3 ! array dimension
real fspl(2,2,inf3,ny) ! spline coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
c
subroutine gridtricub(
> x_newgrid,nx_new,y_newgrid,ny_new,z_newgrid,nz_new,
> f_new,if1,if2,
> nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf4,inf5,iwarn,ier)
c
c regrid a spline function f defined vs. x,y,z as in xpkg, etc.
c to a new grid, given by x_newgrid, y_newgrid, z_newgrid
c
c set warning flag if the range exceeds the range of the
c original x/y/zpkg's.
c
c (xpkg/ypkg/zpkg -- axis data, see genxpkg subroutine)
c
c input:
c
real x_newgrid(nx_new) ! new x grid
real y_newgrid(ny_new) ! new y grid
real z_newgrid(nz_new) ! new z grid
c
c output:
c
integer if1,if2 ! 1st dimensions of f_new
real f_new(if1,if2,nz_new) ! f evaluated on this grid
c
c input:
c
integer nx ! size of old grid
real xpkg(nx,4) ! old grid "package"
integer ny ! size of old grid
real ypkg(ny,4) ! old grid "package"
integer nz ! size of old grid
real zpkg(nz,4) ! old grid "package"
c
integer inf4,inf5 ! array dimensions
real fspl(2,2,2,inf4,inf5,nz) ! spline coefficients of f
c
c output:
c condition codes, =0 for normal exit
c
integer iwarn ! =1 if new grid points out of range
integer ier ! =1 if there is an argument error
c
subroutine fvspline(ict,ivec,ivecd,
> fval,ii,xparam,hx,hxi,fin,nx)
C
C ** use mkspline to set up the fin array **
C
integer ict(3) ! requested output control
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec) ! target cells (i)
real xparam(ivec)
! normalized displacements from (i) corners
C
real hx(ivec) ! grid spacing, and
real hxi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
C
real fin(0:1,nx) ! interpolant data (cf "evspline")
C
real fval(ivecd,3) ! output returned
C
C for detailed description of fin, ict and fval see subroutine
C evspline comments. Note ict is not vectorized; the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii and parameter inputs
C xparam,hx,hxi, are vectorized, and the
C output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
C
subroutine fvbicub(ict,ivec,ivecd,
> fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
> fin,inf2,ny)
C
C ** use mkbicub to set up the fin array **
C
integer ict(6) ! requested output control
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec),jj(ivec) ! target cells (i,j)
real xparam(ivec),yparam(ivec)
! normalized displacements from (i,j) corners
C
real hx(ivec),hy(ivec) ! grid spacing, and
real hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
! & 1/(y(j+1)-y(j))
C
real fin(0:3,inf2,ny) ! interpolant data (cf "evbicub")
C
real fval(ivecd,6) ! output returned
C
C for detailed description of fin, ict and fval see subroutine
C fvbicub comments. Note ict is not vectorized; the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii,jj and parameter inputs
C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
C output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
subroutine fvtricub(ict,ivec,ivecd,
> fval,ii,jj,kk,xparam,yparam,zparam,
> hx,hxi,hy,hyi,hz,hzi,
> fin,inf2,inf3,nz)
C
C use mktricub to set up spline coefficients...
C
integer ict(10) ! requested output control
integer ivec ! vector length
integer ivecd ! vector dimension (1st dim of fval)
C
integer ii(ivec),jj(ivec),kk(ivec) ! target cells (i,j,k)
real xparam(ivec),yparam(ivec),zparam(ivec)
! normalized displacements from (i,j,k) corners
C
real hx(ivec),hy(ivec),hz(ivec) ! grid spacing, and
real hxi(ivec),hyi(ivec),hzi(ivec) ! inverse grid spacing
! 1/(x(i+1)-x(i)) & 1/(y(j+1)-y(j)) & 1/(z(k+1)-z(i))
C
real fin(0:7,inf2,inf3,nz) ! interpolant data (cf "evtricub")
C
real fval(ivecd,10) ! output returned
C
C for detailed description of fin, ict and fval see subroutine
C evtricub comments. Note ict is not vectorized; the same output
C is expected to be returned for all input vector data points.
C
C note that the index inputs ii,jj,kk and parameter inputs
C xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi are vectorized, and
C the output array fval has a vector ** 1st dimension ** whose
C size must be given as a separate argument
C
C to use this routine in scalar mode, pass in ivec=ivecd=1
C
subroutine evspline(xget,x,nx,ilinx,f,ict,fval,ier)
C
C evaluate a 1d cubic Spline interpolant -- this is C2
C
C ** use mkspline to set up the f array **
C
C this subroutine calls two subroutines:
C herm1x -- find cell containing (xget)
C fvspline -- evaluate interpolant function and (optionally)
C derivatives
C
C input arguments:
C ================
C
integer nx ! grid size
real xget ! target of this interpolation
real x(nx) ! ordered x grid
integer ilinx ! ilinx=1 => assume x evenly spaced
C
real f(0:1,nx) ! function data
C
C f(0,i) = f @ x(i)
C f(1,i) = d2f/dx2 @ x(i)
C
C (these are spline coefficients selected for continuous 2-
C diffentiability, see mkspline.for)
C
integer ict(3) ! code specifying output desired
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C ict(3)=1 -- return d2f/dx2 (0, don't)
C
C set ict(1)=3 to get d3f/dx3 (only)
C
C output arguments:
C =================
C
real fval(3) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C fval(3) receives the third output (depends on ict(...) spec)
C
C examples:
C on input ict = [1,1,1]
C on output fval= [f,df/dx,d2f/dx2]
C
C on input ict = [1,0,0]
C on output fval= [f] ... elements 2 -- 3 never referenced.
C
C on input ict = [0,0,1]
C on output fval= [d2f/dx2] ... elements 2 -- 3 never referenced.
C
C ier -- completion code: 0 means OK
subroutine evbicub(xget,yget,x,nx,y,ny,ilinx,iliny,
> f,inf2,ict,fval,ier)
C
C ** use mkbicub to set up the f array **
C
C evaluate a 2d cubic Spline interpolant on a rectilinear
C grid -- this is C2 in both directions.
C
C this subroutine calls two subroutines:
C herm2xy -- find cell containing (xget,yget)
C fvbicub -- evaluate interpolant function and (optionally)
C derivatives
C
C input arguments:
C ================
C
integer nx,ny ! grid sizes
real xget,yget ! target of this interpolation
real x(nx) ! ordered x grid
real y(ny) ! ordered y grid
integer ilinx ! ilinx=1 => assume x evenly spaced
integer iliny ! iliny=1 => assume y evenly spaced
C
real f(0:3,inf2,ny) ! function data
C
C f 2nd dimension inf2 must be .ge. nx
C contents of f:
C
C f(0,i,j) = f @ x(i),y(j)
C f(1,i,j) = d2f/dx2 @ x(i),y(j)
C f(2,i,j) = d2f/dy2 @ x(i),y(j)
C f(3,i,j) = d4f/dx2dy2 @ x(i),y(j)
C
C (these are spline coefficients selected for continuous 2-
C diffentiability, see mkbicub[w].for)
C
integer ict(4) ! code specifying output desired
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C ict(3)=1 -- return df/dy (0, don't)
C ict(4)=1 -- return d2f/dx2 (0, don't)
C ict(5)=1 -- return d2f/dy2 (0, don't)
C ict(6)=1 -- return d2f/dxdy (0, don't)
c the number of non zero values ict(1:6)
c determines the number of outputs...
c
c new dmc December 2005 -- access to higher derivatives (even if not
c continuous-- but can only go up to 3rd derivatives on any one coordinate.
c if ict(1)=3 -- want 3rd derivatives
c ict(2)=1 for d3f/dx3
c ict(3)=1 for d3f/dx2dy
c ict(4)=1 for d3f/dxdy2
c ict(5)=1 for d3f/dy3
c number of non-zero values ict(2:5) gives no. of outputs
c if ict(1)=4 -- want 4th derivatives
c ict(2)=1 for d4f/dx3dy
c ict(3)=1 for d4f/dx2dy2
c ict(4)=1 for d4f/dxdy3
c number of non-zero values ict(2:4) gives no. of outputs
c if ict(1)=5 -- want 5th derivatives
c ict(2)=1 for d5f/dx3dy2
c ict(3)=1 for d5f/dx2dy3
c number of non-zero values ict(2:3) gives no. of outputs
c if ict(1)=6 -- want 6th derivatives
c d6f/dx3dy3 -- one value is returned.
C
C output arguments:
C =================
C
real fval(6) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C fval(3) receives the third output (depends on ict(...) spec)
C fval(4) receives the fourth output (depends on ict(...) spec)
C fval(5) receives the fourth output (depends on ict(...) spec)
C fval(6) receives the fourth output (depends on ict(...) spec)
C
C examples:
C on input ict = [1,1,1,0,0,1]
C on output fval= [f,df/dx,df/dy,d2f/dxdy], elements 5 & 6 not referenced.
C
C on input ict = [1,0,0,0,0,0]
C on output fval= [f] ... elements 2 -- 6 never referenced.
C
C on input ict = [0,0,0,1,1,0]
C on output fval= [d2f/dx2,d2f/dy2] ... elements 3 -- 6 never referenced.
C
C on input ict = [0,0,1,0,0,0]
C on output fval= [df/dy] ... elements 2 -- 6 never referenced.
C
C ier -- completion code: 0 means OK
subroutine evtricub(xget,yget,zget,x,nx,y,ny,z,nz,
> ilinx,iliny,ilinz,
> f,inf2,inf3,ict,fval,ier)
C
C use mktricub to set up spline coefficients...
C
C evaluate a 3d cubic Spline interpolant on a rectilinear
C grid -- this is C2 in all directions.
C
C this subroutine calls two subroutines:
C herm3xyz -- find cell containing (xget,yget,zget)
C fvtricub -- evaluate the spline function (w/derivatives if
C requested).
C
C input arguments:
C ================
C
real xget,yget,zget ! target of this interpolation
real x(nx) ! ordered x grid
real y(ny) ! ordered y grid
real z(nz) ! ordered z grid
integer ilinx ! ilinx=1 => assume x evenly spaced
integer iliny ! iliny=1 => assume y evenly spaced
integer ilinz ! ilinz=1 => assume z evenly spaced
C
real f(0:7,inf2,inf3,nz) ! function data
C
C f 2nd dimension inf2 must be .ge. nx; 3rd dim inf3 .ge. ny
C contents of f:
C
C f(0,i,j,k) = f @ x(i),y(j),z(k)
C f(1,i,j,k) = d2f/dx2 @ x(i),y(j),z(k)
C f(2,i,j,k) = d2f/dy2 @ x(i),y(j),z(k)
C f(3,i,j,k) = d2f/dz2 @ x(i),y(j),z(k)
C f(4,i,j,k) = d4f/dx2dy2 @ x(i),y(j),z(k)
C f(5,i,j,k) = d4f/dx2dz2 @ x(i),y(j),z(k)
C f(6,i,j,k) = d4f/dy2dz2 @ x(i),y(j),z(k)
C f(7,i,j,k) = d6f/dx2dy2dz2 @ x(i),y(j),z(k)
C
integer ict(10) ! code specifying output desired
C
C ict(1)=1 -- return f (0, don't)
C ict(2)=1 -- return df/dx (0, don't)
C ict(3)=1 -- return df/dy (0, don't)
C ict(4)=1 -- return df/dz (0, don't)
C ict(5)=1 -- return d2f/dx2 (0, don't)
C ict(6)=1 -- return d2f/dy2 (0, don't)
C ict(7)=1 -- return d2f/dz2 (0, don't)
C ict(8)=1 -- return d2f/dxdy (0, don't)
C ict(9)=1 -- return d2f/dxdz (0, don't)
C ict(10)=1-- return d2f/dydz (0, don't)
c
c (new dmc Dec 2005 -- higher derivatives available)
c ict(1)=3 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:8) select: fxxy fxxz fxyy fxyz fxzz fyyz fyzz
c ->note ict(1)=3, ict(5)=1 gives fxyz = d3f/dxdydz
c ict(1)=-3 --> 3rd derivative, 3 in one coordinate
c ict(2:4) select: fxxx fyyy fzzz
c ict(1)=4 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:7) select: fxxyy fxxyz fxxzz fxyyz fxyzz fyyzz
c ict(1)=-4 --> 3rd derivative, 3 in one coordinate
c ict(2:7) select: fxxxy fxxxz fxyyy fxzzz fyyyz fyzzz
c ict(1)=5 --> 3rd derivative, .le.2 diff. in any coordinate
c ict(2:4) select: fxxyyz fxxyzz fxyyzz
c ict(1)=-5 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyy fxxxyz fxxxzz fxxyyy fxxzzz
c fxyyyz fxyzzz fyyyzz fzzzyy
c ict(1)=6 --> 3rd derivative, .le.2 diff. in any coordinate
c fxxyyzz
c ict(1)=-6 --> 3rd derivative, 3 in one coordinate
c ict(2:10) select: fxxxyyy fxxxyyz fxxxyzz fxxxyyz
c fxxyyyz fxxyzzz fxyyyzz fxyyzzz fyyyzzz
c ict(1)=-7 --> 7th derivative
c ict(2:7) select: fxxxyyyz fxxxyyzz fxxxyzzz
c fxxyyyzz fxxyyzzz fxyyyzzz
c ict(1)=-8 --> 8th derivative
c ict(2:4) select: fxxxyyyzz fxxxyyzzz fxxyyyzzz
c ict(1)=-9 --> 9th derivative: fxxxyyyzzz
C
C output arguments:
C =================
C
real fval(10) ! output data
integer ier ! error code =0 ==> no error
C
C fval(1) receives the first output (depends on ict(...) spec)
C fval(2) receives the second output (depends on ict(...) spec)
C fval(3) receives the third output (depends on ict(...) spec)
C fval(4) receives the 4th output (depends on ict(...) spec)
C fval(5-10) receive 5th thru 10th outputs (if required by ict(...) spec)
C
C examples:
C on input ict = [1,1,1,1,0,0,0,0,0,0,0]
C on output fval= [f,df/dx,df/dy,df/dz]
C
C on input ict = [1,0,0,0,0,0,0,0,0,0,0]
C on output fval= [f] ... elements 2-10 never referenced
C
C on input ict = [0,1,1,0,0,0,0,0,0,0,0]
C on output fval= [df/dx,df/dy] ... elements 3-10 never referenced
C
C on input ict = [0,0,0,0,1,0,0,0,0,0,0]
C on output fval= [d2f/dx2] ... elements 2-10 never referenced.
C
C ier -- completion code: 0 means OK
PSPLINE's "explicit splines" are somewhat faster to evaluate than
PSPLINE's "compact splines", but they require far more memory.
Spline setup routines:
explicit spline routines memory requirement
------------------------------------------------------------
1d: cspline 4*Nx
2d: bcspline 16*Nx*Ny
3d: tcspline 64*Nx*Ny*Nz
==> these routines create the spline coefficients and any other
information needed by the lookup / evaluation routines --
with various boundary condition control options.
Vectorized Evaluation Routines (recommended):
=============================================
given a location vector, do lookup & evaluation:
1d 2d 3d
spvec bcspvec tcspvec
given a grid (different than the original grid), to
map the function from the old grid to the new grid:
spgrid bcspgrid tcspgrid
if zone location data is already computed, to do just
the vectorized spline evaluation:
cspevfn bcspevfn tcspevfn
Scalar Routines:
================
Routines for compact spline function evaluation (scalar, grid
lookup included):
cspeval bcspeval tcspeval
Generally, the spline setup routines require workspaces. As these
routines are written in highly portable fortran-77 code (and thus
lack dynamic memory allocation capability), calling programs need
to declare and provide writable workspace arrays.
The following table summarizes the sizes of workspaces needed.
In this table, "x" denotes the 1st dimension, "y" the 2nd
dimension, and "z" the 3rd dimension of the input data.
=> 1d spline routines:
cspline workspace size: 1*nx (periodic spline option only; for
other options the workspace is unused and a dummy
argument can be provided).
=> 2d spline routines:
bcspline: ws size: 5*max(nx,ny) --or--
5*max(nx,ny)+4*nx*ny [note#1]
=> 3d spline routines:
tcspline: ws size: 10*max(nx,ny,nz) + 20*nx*ny --or--
16*nx*ny*nz [note#2]
[note#1] the larger number is required if the "y" boundary
condition set is "non-linear" (see below).
[note#2] the larger number is required if the "z" boundary
condition set is "non-linear".
A boundary condition set BC is considered "linear" if and only
if, given two splines s1 and s2 satisfying BC, the spline given
by the formula
s = a*s1 + s2 (a is a real number)
also satisfies BC.
Examples of linear boundary condition sets:
* periodic: s'(x(1))=s'(x(nx)) & s''(x(1))=s''(x(nx))
or any combination of two of the following:
* not a knot: s'''(x) continuous at x(2) or x(nx-1).
* zero derivative: s'(x(1))=0 and/or s'(x(nx))=0
* zero 2nd derivative: s''(x(1))=0 and/or s''(x(nx))=0
* a derivative estimated by a divided difference formula
Non-linear boundary condition sets are those for which an
explicit non-zero derivative or 2nd derivative value is
specified. Splines satisfying such boundary conditions
cannot be added without alteration of the boundary condition.
c cspline -- dmc 15 Feb 1999
c
c a standard interface to various 1d spline setup routines
c
subroutine cspline(x,nx,fspl,ibcxmin,bcxmin,ibcxmax,bcxmax,
> wk,iwk,ilinx,ier)
c
real x(nx) ! x axis (in)
real fspl(4,nx) ! spline data (in/out)
integer ibcxmin ! x(1) BC flag (in, see comments)
real bcxmin ! x(1) BC data (in, see comments)
integer ibcxmax ! x(nx) BC flag (in, see comments)
real bcxmax ! x(nx) BC data (in, see comments)
real wk(iwk) ! workspace of size at least nx
c
c workspace only used if "periodic" boundary condition is selected.
c
integer ilinx ! even spacing flag (out)
integer ier ! output status, =0 means OK
c
c this routine computes spline coefficients for a 1d spline --
c evaluation of the spline can be done by cspeval.for subroutines
c or directly by inline code.
c
c the input x axis x(1...nx) must be strictly ascending, i.e.
c x(i+1).gt.x(i) is required for i=1 to nx-1. This is checked and
c ier=1 is set and the routine exits if the test is not satisfied.
c
c on output, ilinx=1 is set if, to a reasonably close tolerance,
c all grid spacings x(i+1)-x(i) are equal. This allows a speedier
c grid lookup algorithm on evaluation of the spline. If on output
c ilinx=2, this means the spline x axis is not evenly spaced.
c
c the input data for the spline are given in f[j] = fspl(1,j). The
c output data are the spline coefficients fspl(2,j),fspl(3,j), and
c fspl(4,j), j=1 to nx. The result is a spline s(x) satisfying the
c boundary conditions and with the properties
c
c s(x(j)) = f(1,j)
c s'(x) is continuous even at the grid points x(j)
c s''(x) is continuous even at the grid points x(j)
c
c the formula for evaluation of s(x) is:
c
c let dx = x-x(i), where x(i).le.x.le.x(i+1). Then,
c s(x)=f(1,i) + dx*(f(2,i) +dx*(f(3,i) + dx*f(4,i)))
c
c ==>boundary conditions. Complete specification of a 1d spline
c requires specification of boundary conditions at x(1) and x(nx).
c
c this routine provides 4 options:
c
c ibcxmin=-1 -- periodic boundary condition. This means the
c boundary conditions s'(x(1))=s'(x(nx)) and s''(x(1))=s''(x(nx))
c are imposed. Note that s(x(1))=s(x(nx)) (i.e. f(1,1)=f(1,nx))
c is not required -- that is determined by the f array input data.
c The periodic boundary condition is to be preferred for periodic
c data. When splining periodic data f(x) with period P, the relation
c x(nx)=x(1)+n*P, n = the number of periods (usually 1), should hold.
c (ibcxmax, bcxmin, bcxmax are ignored).
c
c if a periodic boundary condition is set, this covers both boundaries.
c for the other types of boundary conditions, the type of condition
c chosen for the x(1) boundary need not be the same as the type chosen
c for the x(nx) boundary.
c
c ibcxmin=0 | ibcxmax=0 -- this specifies a "not a knot" boundary
c condition -- see cubsplb.for. This is a common way for inferring
c a "good" spline boundary condition automatically from data in the
c vicinity of the boundary. (bcxmin | bcxmax are ignored).
c
c ibcxmin=1 | ibcxmax=1 -- boundary condition is to have s'(x(1)) |
c s'(x(nx)) match the passed value (bcxmin | bcxmax).
c
c ibcxmin=2 | ibcxmax=2 -- boundary condition is to have s''(x(1)) |
c s''(x(nx)) match the passed value (bcxmin | bcxmax).
c
c ibcxmin=3 | ibcxmax=3 -- boundary condition is to have s'(x(1)) |
c s'(x(nx)) = 0.0 (ZERO)
c
c ibcxmin=4 | ibcxmax=4 -- boundary condition is to have s''(x(1)) |
c s''(x(nx)) = 0.0 (ZERO)
c
c ibcxmin=5 | ibcxmax=5 == BC is to have s' match the 1st divided
c difference at the endpoint.
c
c ibcxmin=6 | ibcxmax=6 == BC is to have s'' match the 2nd divided
c difference at the endpoint.
c
c ibcxmin=7 | ibcxmax=7 == BC is to have s''' match the 3rd divided
c difference at the endpoint.
c
c set up coefficients for bicubic spline with following BC's:
c FULL BC CONTROL at all bdys
c
subroutine bcspline(x,inx,th,inth,fspl,inf3,
> ibcxmin,bcxmin,ibcxmax,bcxmax,
> ibcthmin,bcthmin,ibcthmax,bcthmax,
> wk,nwk,ilinx,ilinth,ier)
c
real x(inx),th(inth),fspl(4,4,inf3,inth),wk(nwk)
real bcxmin(inth),bcxmax(inth)
real bcthmin(inx),bcthmax(inx)
c
c input:
c x(1...inx) -- abscissae, first dimension of data
c th(1...inth) -- abscissae, second (periodic) dimension of data
c fspl(1,1,1..inx,1..inth) -- function values
c inf3 -- fspl dimensioning, inf3.ge.inx required.
c
c boundary conditions input:
c ibcxmin -- indicator for boundary condition at x(1):
c bcxmin(...) -- boundary condition data
c =-1 -- periodic boundary condition
c =0 -- use "not a knot", bcxmin(...) ignored
c =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
c =2 -- match 2nd derivative, specified at x(1),th(ith) by
c bcxmin(ith)
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
c =5 -- match 1st derivative to 1st divided difference
c =6 -- match 2nd derivative to 2nd divided difference
c =7 -- match 3rd derivative to 3rd divided difference
c (for more detailed definition of BCs 5-7, see the
c comments of subroutine mkspline)
c the bcxmin(...) array is referenced only if ibcxmin=1 or
c ibcxmin=2.
c
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c NOTE: if ibcxmin=-1, ibcxmax is ignored! ...and the BC is periodic.
c
c ibcthmin -- indicator for boundary condition at th(1):
c bcthmin(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c ibcthmax -- indicator for boundary condition at th(inth):
c bcthmax(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c NOTE: if ibcthmin=-1, ibcthmax is ignored! ...the BC is periodic.
c
c NOTE the bcxmin,bcxmax,bcthmin,bcthmax arrays are only used if the
c corresponding boundary condition flags are set to 1 or 2.
c Carefully note the dimensioning of these arrays!
c
c output:
c fspl(*,*,1..inx,1..inth) -- bicubic spline coeffs (4x4)
c ...fspl(1,1,*,*) is not replaced.
c
c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced
c (tol=1e-3)
c ilinth-- =1 on output if th(inth) are nearly evenly spaced
c (tol=1e-3)
c
c ier -- completion code, 0 for normal
c
c workspace:
c wk -- must be at least 5*max(inx,inth) large *** OR ***
c 4*inx*inth + 5*max(inx,inth) if ibcthmin=1 or 2
c or ibcthmax=1 or 2 and non-zero explicit derivative
c boundary information is provided.
c nwk -- size of workspace provided.
c
c tcspline -- dmc 20 Jan 1999
c
c set up coefficients for bicubic spline with following BC's:
c * LHS and RHS handled as in cubspl.for for 1st coordinate
c * derivatives periodic in second coordinate (use pspline.for)
c
subroutine tcspline(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
> ibcxmin,bcxmin,ibcxmax,bcxmax,inb1x,
> ibcthmin,bcthmin,ibcthmax,bcthmax,inb1th,
> ibcphmin,bcphmin,ibcphmax,bcphmax,inb1ph,
> wk,nwk,ilinx,ilinth,ilinph,ier)
c
real x(inx),th(inth),ph(inph)
real fspl(4,4,4,inf4,inf5,inph),wk(nwk)
real bcxmin(inb1x,inph),bcxmax(inb1x,inph) ! inth x inph used
real bcthmin(inb1th,inph),bcthmax(inb1th,inph) ! inx x inph used
real bcphmin(inb1ph,inth),bcphmax(inb1ph,inth) ! inx x inth used
c
c input:
c x(1...inx) -- abscissae, first dimension of data
c th(1...inth) -- abscissae, second (periodic) dimension of data
c ph(1...inph) -- abscissae, third (periodic) dimension of data
c fspl(1,1,1,1..inx,1..inth,1..inph) -- function values
c inf4 -- fspl dimensioning, inf4.ge.inx required.
c inf5 -- fspl dimensioning, inf5.ge.inth required.
c
c boundary conditions input:
c
c bc data at xmin, xmax vs. theta,phi
c bc data at thmin, thmax vs. x,phi
c bc data at phmin, phmax vs. x,theta
c
c ibcxmin -- indicator for boundary condition at x(1):
c bcxmin(...) -- boundary condition data
c =-1 -- use periodic boundary condition
c =0 -- use "not a knot", bcxmin(...) ignored
c =1 -- match slope, specified at x(1),th(ith),ph(iph)
c by bcxmin(ith,iph)
c =2 -- match 2nd derivative, specified at x(1),th(ith),ph(iph)
c by bcxmin(ith,iph)
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all pts
c =4 -- boundary condition is d2f/dx2=0 at x(1), all (th,ph) pts
c =5 -- match 1st derivative to 1st divided difference
c =6 -- match 2nd derivative to 2nd divided difference
c =7 -- match 3rd derivative to 3rd divided difference
c (for more detailed definition of BCs 5-7, see the
c comments of subroutine mkspline)
c bcxmin(...) is referenced ONLY if ibxcmin=1 or ibxcmin=2.
c
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpretation as with ibcxmin, bcxmin)
c
c NOTE: if ibcxmin=-1 then the periodic BC applies on both sides
c and ibcxmax is ignored.
c inb1x -- 1st dimension of bcxmin, bcxmax: if ibcxmin .gt. 0
c or ibcxmax .gt. 0, this must be .ge. inth.
c
c interpretation of ibcthmin,bcthmin,ibcthmax,bcthmax,inb1th
c is same as with ibcxmin,...
c
c interpretation of ibcphmin,bcphmin,ibcphmax,bcphmax,inb1ph
c is same as with ibcxmin,...
c
c output:
c fspl(*,*,*,1..inx,1..inth,1..inph) -- bicubic spline coeffs (4x4)
c ...fspl(1,1,1,*,*,*) is not replaced.
c
c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced
c (tol=1e-3)
c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
c ilinph-- =1 on output if ph(inph) evenly spaced (tol=1e-3)
c
c ier -- completion code, 0 for normal
c
c workspace:
c wk -- because tcspline uses bcspline as a subroutine, a large
c workspace is needed, depending on boundary conditions:
c
c if ibcphmin = 1 or 2, or, ibcphmax = 1 or 2, and non-zero
c explicit d/dphi or d2/dphi2 boundary condition data is
c provided, the workspace size must be at least
c
c 16*inx*inth*inph.
c
c if ibcphmin and ibcphmax .ne. 1 or 2, and no explicit
c non-zero d/dphi or d2/dphi2 boundary condition data is
c provided, then, the workspace size must be at least:
c
c 20*inx*inth + 10*max(inx,inth,inph).
c
c nwk -- actual size of workspace
c
There is a small collection of routines which have fewer arguments,
but also less control over boundary conditions. These are shown here.
c
c set up coefficients for bicubic spline with following BC's:
c * LHS and RHS handled as in cubspl.for for 1st coordinate
c * derivatives periodic in second coordinate (use pspline.for)
c
subroutine bpsplinb(x,inx,th,inth,fspl,inf3,
> ibcxmin,bcxmin,ibcxmax,bcxmax,
> wk,nwk,ilinx,ilinth,ier)
c
real x(inx),th(inth),fspl(4,4,inf3,inth),wk(nwk)
real bcxmin(inth),bcxmax(inth)
c
c input:
c x(1...inx) -- abscissae, first dimension of data
c th(1...inth) -- abscissae, second (periodic) dimension of data
c fspl(1,1,1..inx,1..inth) -- function values
c inf3 -- fspl dimensioning, inf3.ge.inx required.
c
c boundary conditions input:
c ibcxmin -- indicator for boundary condition at x(1):
c bcxmin(...) -- boundary condition data
c =0 -- use "not a knot", bcxmin(...) ignored
c =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
c =2 -- match 2nd derivative, specified at x(1),th(ith) by
c bcxmin(ith)
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
c =5 -- match 1st derivative to 1st divided difference
c =6 -- match 2nd derivative to 2nd divided difference
c =7 -- match 3rd derivative to 3rd divided difference
c (for more detailed definition of BCs 5-7, see the
c comments of subroutine mkspline)
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpolation as with ibcxmin, bcxmin)
c
c output:
c fspl(*,*,1..inx,1..inth) -- bicubic spline coeffs (4x4)
c ...fspl(1,1,*,*) is not replaced.
c
c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced
c (tol=1e-3)
c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
c
c ier -- completion code, 0 for normal
c
c workspace:
c wk -- must be at least 5*max(inx,inth) large
c nwk -- size of workspace
c
c set up coefficients for bicubic spline with following BC's:
c * LHS and RHS handled as in spline.for for 1st coordinate
c * derivatives periodic in second coordinate (use pspline.for)
c
c see similar routine, bpsplinb, for more control over x boundary
c conditions...
c
subroutine bpspline(x,inx,th,inth,fspl,inf3,
> wk,nwk,ilinx,ilinth,ier)
c
real x(inx),th(inth),fspl(4,4,inf3,inth),wk(nwk)
integer ilinx,ilinth
c
c input:
c x(1...inx) -- abscissae, first dimension of data
c th(1...inth) -- abscissae, second (periodic) dimension of data
c fspl(1,1,1..inx,1..inth) -- function values
c inf3 -- fspl dimensioning, inf3.ge.inx required.
c
c output:
c fspl(*,*,1..inx,1..inth) -- bicubic spline coeffs (4x4)
c ...fspl(1,1,*,*) is not replaced.
c
c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced
c (tol=1e-3)
c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
c
c ier -- completion code, 0 for normal
c
c workspace:
c wk -- must be at least 5*max(inx,inth) large
c nwk -- size of workspace
c
c set up coefficients for bicubic spline with following BC's:
c * LHS and RHS handled as in cubspl.for for 1st coordinate
c * derivatives periodic in second coordinate (use pspline.for)
c
subroutine tpsplinb(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
> ibcxmin,bcxmin,ibcxmax,bcxmax,inb1,
> wk,nwk,ilinx,ilinth,ilinph,ier)
c
real x(inx),th(inth),ph(inph)
real fspl(4,4,4,inf4,inf5,inph),wk(nwk)
real bcxmin(inb1,inph),bcxmax(inb1,inph)
c
c input:
c x(1...inx) -- abscissae, first dimension of data
c th(1...inth) -- abscissae, second (periodic) dimension of data
c ph(1...inph) -- abscissae, third (periodic) dimension of data
c fspl(1,1,1,1..inx,1..inth,1..inph) -- function values
c inf4 -- fspl dimensioning, inf4.ge.inx required.
c inf5 -- fspl dimensioning, inf5.ge.inth required.
c
c boundary conditions input:
c ibcxmin -- indicator for boundary condition at x(1):
c bcxmin(...) -- boundary condition data
c =0 -- use "not a knot", bcxmin(...) ignored
c =1 -- match slope, specified at x(1),th(ith),ph(iph)
c by bcxmin(ith,iph)
c =2 -- match 2nd derivative, specified at x(1),th(ith),ph(iph)
c by bcxmin(ith,iph)
c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
c =5 -- match 1st derivative to 1st divided difference
c =6 -- match 2nd derivative to 2nd divided difference
c =7 -- match 3rd derivative to 3rd divided difference
c (for more detailed definition of BCs 5-7, see the
c comments of subroutine mkspline)
c ibcxmax -- indicator for boundary condition at x(nx):
c bcxmax(...) -- boundary condition data
c (interpolation as with ibcxmin, bcxmin)
c
c output:
c fspl(*,*,*,1..inx,1..inth,1..inph) -- bicubic spline coeffs (4x4)
c ...fspl(1,1,1,*,*,*) is not replaced.
c
c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced
c (tol=1e-3)
c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
c ilinph-- =1 on output if ph(inph) evenly spaced (tol=1e-3)
c
c ier -- completion code, 0 for normal
c
c workspace:
c wk -- must be at least 5*max(inx,inth,inph) large
c nwk -- size of workspace
c set up coefficients for bicubic spline with following BC's:
c * 1st dimension: bounding cubic (divided difference)
c (use spline.for)
c * derivatives periodic in second coordinate (use pspline.for)
c * derivatives periodic in third coordinate (use pspline.for)
c
c for more control over boundary conditions, use tpspline or tcspline.
c
c for evaluation of interpolant, see tcspeval.for
c
subroutine tpspline(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
> wk,nwk,ilinx,ilinth,ilinph,ier)
c
real x(inx),th(inth),ph(inph)
real fspl(4,4,4,inf4,inf5,inph),wk(nwk)
c
c input:
c x(1...inx) -- abscissae, first dimension of data
c th(1...inth) -- abscissae, second (periodic) dimension of data
c ph(1...inph) -- abscissae, third (periodic) dimension of data
c fspl(1,1,1,1..inx,1..inth,1..inph) -- function values
c inf4 -- fspl dimensioning, inf4.ge.inx required.
c inf5 -- fspl dimensioning, inf5.ge.inth required.
c
c output:
c fspl(*,*,*,1..inx,1..inth,1..inph) -- bicubic spline coeffs (4x4)
c ...fspl(1,1,1,*,*,*) is not replaced.
c
c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced
c (tol=1e-3)
c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
c ilinph-- =1 on output if ph(inph) evenly spaced (tol=1e-3)
c
c ier -- completion code, 0 for normal
c
c workspace:
c wk -- must be at least 5*max(inx,inth,inph) large
c nwk -- size of workspace
c
subroutine spvec(ict,ivec,xvec,ivd,fval,nx,xpkg,fspl,iwarn,ier)
c
c vectorized spline evaluation routine -- 1d spline
c 1. call vectorized zone lookup routine
c 2. call vectorized spline evaluation routine
c
c--------------------------
c input:
integer ict(3) ! selector:
c ict(1)=1 for f (don't evaluate if ict(1)=0)
c ict(2)=1 for df/dx (don't evaluate if ict(2)=0)
c ict(3)=1 for d2f/dx2 (don't evaluate if ict(3)=0)
c
c set ict(1)=3 to get d3f/dx3 (only)
c
integer ivec ! vector dimensioning
c
c ivec-- number of vector pts (spline values to look up)
c
real xvec(ivec) ! x-locations at which to evaluate
c
integer ivd ! 1st dimension of output array
c
c ivd -- 1st dimension of fval, .ge.ivec
c
c output:
real fval(ivd,*) ! output array
c
c fval(1:ivec,1) -- values