PSPLINE

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.

Home | Top |

Feb_2010_Improvement

Many routines in pspline used to insist on a minimum gridsize of nx=4. It is now allowed to use nx=2 (the smallest conceivable size, containing only a single cubic or linear element). The test codes (e.g. pspltsub in pspltest) have been extended with a new subroutine "smalltest" that verifies handling of boundary condition controls in spline routines with small grids (some special case coding was required for implementation).

Home | Top |

Feb_2004_Fixes_and_Improvements

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.

Home | Top |

Acknowledgement_and_details

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.

Home | Top |

Acknowledgements

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.

Home | Top |

Case_study

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.

Home | Top |

Spline_Basics

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

Home | Top |

Bicubic_Splines

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.

Home | Top |

Existence_and_Uniqueness

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.

Home | Top |

Construction

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.

Home | Top |

Compact_Algorithm

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.

Home | Top |

NAG_comparison

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.

Home | Top |

Tricubic_Splines

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).

Home | Top |

Higher_dimensionalities

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.

Home | Top |

Hermite_interpolation_basics

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

Home | Top |

Bicubic_Hermite

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.

Home | Top |

Hermite_setup_2d

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.

Home | Top |

Tricubic_Hermite

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.

Home | Top |

Hermite_setup_3d

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.

Home | Top |

Performance_considerations

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)!

Home | Top |

Quality_of_results

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]

Home | Top |

Recommended_strategy

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.

Home | Top |

EZspline_software_interface

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

Home | Top |

F77_standard_software

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.

Home | Top |

interface_module

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.

Home | Top |

genxpkg

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.

Home | Top |

summary_of_lookup_test_results

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)

Home | Top |

brief_description_of_algorithms

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.

Home | Top |

xlookup

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!

Home | Top |

Compact_Splines

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

Home | Top |

mkspline

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.

Home | Top |

mkbicub

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

Home | Top |

mktricub

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

Home | Top |

vecspline

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

Home | Top |

vecbicub

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

Home | Top |

vectricub

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

Home | Top |

gridspline

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

Home | Top |

gridbicub

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

Home | Top |

gridtricub

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

Home | Top |

fvspline

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

Home | Top |

fvbicub

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

Home | Top |

fvtricub

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

Home | Top |

evspline

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

Home | Top |

evbicub

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

Home | Top |

evtricub

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

Home | Top |

Explicit_Splines

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

Home | Top |

workspace_arrays

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.

Home | Top |

cspline

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

Home | Top |

bcspline

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

Home | Top |

tcspline

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

Home | Top |

legacy_setup_routines

There is a small collection of routines which have fewer arguments, but also less control over boundary conditions. These are shown here.

Home | Top |

bpsplinb

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

Home | Top |

bpspline

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

Home | Top |

tpsplinb

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

Home | Top |

tpspline

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

Home | Top |

spvec

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 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(4,nx) ! (non-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

Home | Top |

bcspvec

subroutine bcspvec(ict,ivec,xvec,yvec,ivd,fval, > nx,xpkg,ny,ypkg,fspl,inf3, > iwarn,ier) c c vectorized spline evaluation routine -- 2d 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 inf3 ! fspl 3rd array dimension, .ge.nx real fspl(4,4,inf3,ny) ! (non-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

Home | Top |

tcspvec

subroutine tcspvec(ict,ivec,xvec,yvec,zvec,ivd,fval, > nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf4,inf5, > iwarn,ier) c c vectorized spline evaluation routine -- 3d 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 inf4 ! fspl 4th array dimension, .ge.nx integer inf5 ! fspl 5th array dimension, .ge.ny real fspl(4,4,4,inf4,inf5,nz) ! (non-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

Home | Top |

spgrid

subroutine spgrid(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(4,nx) ! 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

Home | Top |

bcspgrid

subroutine bcspgrid( > 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 of x_newgrid or y_newgrid exceeds c the range of the original xpkg or ypkg. c c (xpkg,ypkg arrays -- 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(4,4,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

Home | Top |

tcspgrid

subroutine tcspgrid( > 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,ypkg,zpkg c to a new grid, given by x_newgrid, y_newgrid, z_newgrid c c set warning flag if the range of x_newgrid or y_newgrid or z_newgrid c exceeds the range of the original xpkg or ypkg or zpkg. c c (xpkg,ypkg,zpkg arrays -- 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 ! array dimension integer inf5 ! array dimension real fspl(4,4,4,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

Home | Top |

cspevfn

==> note, in real life, if speed is a concern, the caller should evaluate the spline directly, using the coefficients array. c cspevfn -- OK now evaluate the cubic spline c subroutine cspevfn(ict,ivec,ivd,fval,iv,dxv,f,nx) c c input: integer ict(3) ! selector: c ict(1)=1 for f (don't evaluate f if ict(1)=0) c ict(2)=1 for df/dx "" c ict(3)=1 for d2f/dx2 c c set ict(1)=3 to get d3f/dx3 (only) c integer ivec,ivd ! vector dimensioning c c ivec-- number of vector pts (spline values to look up) c ivd -- 1st dimension of fval, .ge.ivec c c output: real fval(ivd,6) ! output array c c v = index to element in vector; c fval(v,1) = first item requested by ict(...), c fval(v,2) = 2nd item requested, ...etc... c c input: integer iv(ivec) ! grid cell indices -- vectors real dxv(ivec) ! displacements w/in cell -- vectors c real f(4,nx) ! cubic fcn spline coeffs array c c usage example: c c 1. for each element xx(v) in a vector of x values: c find the x zone index and displacement with respect to the c lower end point of the zone; store thes in vectors iv and dxv. c c 2. set ict(1)=0, ict(2)=1, ict(3)=0 -- get only the 1st derivative c c 3. ivec is the length of the vector; ivd is the 1st dimension of c the array fval to receive the output c c real fval(ivd,3) c real xv(ivd) c real dxv(ivd) c integer iv(ivd) c integer ict(3) c c real fspline(4,nx) ! spline coeffs c data ict/0,1,0/ ! this call: want 1st derivative only c ! this will be output to fval(*,1). c ... c do iv=1,ivec c ... ! compute indices and displacements c enddo c call cspevfn(ict,ivec,ivd,fval,iv,dxv,fspline,nx) c

Home | Top |

bcspevfn

c bcspevfn -- OK now evaluate the bicubic spline c subroutine bcspevfn(ict,ivec,ivd,fval,iv,jv,dxv,dyv,f,inf3,ny) c c input: integer ict(6) ! selector: c ict(1)=1 for f (don't evaluate f if ict(1)=0) c ict(2)=1 for df/dx "" c ict(3)=1 for df/dy "" c ict(4)=1 for d2f/dx2 c ict(5)=1 for d2f/dy2 c ict(6)=1 for d2f/dxdy c 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 integer ivec,ivd ! vector dimensioning c c ivec-- number of vector pts (spline values to look up) c ivd -- 1st dimension of fval, .ge.ivec c c output: real fval(ivd,6) ! output array c c v = index to element in vector; c fval(v,1) = first item requested by ict(...), c fval(v,2) = 2nd item requested, ...etc... c c input: integer iv(ivec),jv(ivec) ! grid cell indices -- vectors real dxv(ivec),dyv(ivec) ! displacements w/in cell -- vectors c real f(4,4,inf3,ny) ! bicubic fcn spline coeffs array integer inf3 ! 3rd dimension of f -- .ge. nx c c usage example: c c 1. for each element (xx(v),yy(v)) in a vector of (x,y) pairs, c find the x and y zone indices and displacements with respect c to the "lower left corner" of the zone; store these in vectors c iv,jv and dxv,dyv. c c 2. set ict(1)=0, ict(2)=1, ict(3)=1, the rest zero -- get only c the 1st derivatives. c c 3. ivec is the length of the vector; ivd is the 1st dimension c of the array fval to receive the output c c real fval(ivd,6) c real xv(ivd),yv(ivd) c integer iv(ivd),jv(ivd) c real dxv(ivd),dyv(ivd) c integer ict(6) c c real fspline(4,4,nx,ny) ! spline coeffs c data ict/0,1,1,0,0,0/ ! this call: want 1st derivatives c ! only ... these will be output to c ! fval(*,1) fval(*,2) c ... c do iv=1,ivec c ... ! find indices and displacements c enddo c call bcspevfn(ict,ivec,ivd,fval,iv,jv,dxv,dyv,fspline,nx,ny) c

Home | Top |

tcspevfn

c tcspevfn -- OK now evaluate the tricubic spline c evaluate at a vector set of target locations as specified by c input vectors (iv,jv,kv), (dxv,dyv,dzv) c subroutine tcspevfn(ict,ivec,ivd,fval,iv,jv,kv,dxv,dyv,dzv, > f,inf4,inf5,nz) c c input: integer ict(10) ! selector: c ict(1)=1 for f (don't evaluate f if ict(1)=0) c ict(2)=1 for df/dx "" c ict(3)=1 for df/dy "" c ict(4)=1 for df/dz "" c ict(5)=1 for d2f/dx2 c ict(6)=1 for d2f/dy2 c ict(7)=1 for d2f/dz2 c ict(8)=1 for d2f/dxdy c ict(9)=1 for d2f/dxdz c ict(10)=1 for d2f/dydz 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,ivd ! vector dimensioning c c ivec-- number of vector pts (spline values to look up) c ivd -- 1st dimension of fval, .ge. ivec c c output: real fval(ivd,10) ! output array c c for vector elements v, (iv(v),jv(v),kv(v),dxv(v),dyv(v),dzv(v)) c c fval(v,1) = first item requested by ict(...), c fval(v,2) = 2nd item requested, ...etc... c c input: integer iv(ivec),jv(ivec),kv(ivec) ! grid cell indices real dxv(ivec),dyv(ivec),dzv(ivec) ! displacements w/in cell c real f(4,4,4,inf4,inf5,nz) ! tricubic fcn spline coeffs array integer inf4 ! 4th dimension of f, .le.nx integer inf5 ! 5th dimension of f, .le.ny c c usage example: c c 1. for each element (xx(v),yy(v),zz(v)) in a vector of (x,y,z) c triples, find the x,y,z zone indices and displacements with c to the "lower left corner" of the zone; store these in vectors c iv,jv,kv and dxv,dyv,dzv c c 2. set ict(1)=0, ict(2)=1, ict(3)=1, ict(4)=1 & the rest zero -- c to get only the 1st derivatives. c c 3. ivec is the length of the vector; ivd is the 1st dimension c of the array fval to receive the output c c real fval(ivd,10) c real xv(ivd),yv(ivd),zv(ivd) c integer iv(ivd),jv(ivd),kv(ivd) c real dxv(ivd),dyv(ivd),dzv(ivd) c integer ict(10) c c real fspline(4,4,4,nx,ny,nz) ! spline coeffs c data ict/0,1,1,1,0,0,0,0,0,0/ ! this call: want 1st c ! derivatives only ... these will c ! be output to c ! fval(*,1) fval(*,2) fval(*,3) c ... c do iv=1,ivec c ... ! find indices and displacements c enddo c call tcspevfn(ict,ivec,ivd,fval,iv,jv,kv,dxv,dyv,dzv, c > fspline,nx,ny,nz) c

Home | Top |

cspeval

c cspeval -- eval cubic spline function and/or derivatives c subroutine cspeval(xget,iselect,fval,x,nx,ilinx,f,ier) c real xget ! interpolation target real fval(3) ! output values real x(nx),f(4,nx) ! spline data c integer iselect(3) ! output selector integer ilinx ! =1 if x(...) is evenly spaced c c input: c (xget) location where interpolated value is desired c x(1).le.xget.le.x(nx) expected c c iselect select desired output c c iselect(1)=1 -- want function value (f) itself c iselect(2)=1 -- want df/dx c iselect(3)=1 -- want d2f/dx2 c c example: iselect(1)=iselect(2)=iselect(3)=1 c f, df/dx, and d2f/dx2 are all evaluated c and returned in order in fval(1), fval(2), c and fval(3) c iselect(1)=0, iselect(2)=1, iselect(3)=0 c only the 1st derivative is evaluated c and returned in fval(1). c c set iselect(1)=3 to get d3f/dx3, 1 value only. c c see fval (output) description. c c x(1...nx) independent coordinate x, strict ascending c c ilinx -- =1: flag that x is linearly spaced c (find zone at once, avoid search for speed) c c **CAUTION** actual even spacing of x, is NOT CHECKED HERE! c c f the function values (at grid points) and spline c coefficients c c evaluation formula: for point x btw x(i) and x(i+1), dx=x-x(i) c c spline value = c f(1,i) + dx*f(2,i) + dx**2*f(3,i) + dx**3*f(4,i) c c output: c up to 3 elements of fval, ordered as follows: c fval(1)=function value or lowest order derivative requested c fval(2)=next order derivative c etc c the ordering is a subset of the sequence given under the c "iselect" description. c c ier = 0 -- successful completion; = 1 -- an error occurred.

Home | Top |

bcspeval

c bcspeval -- eval bicubic spline function and/or derivatives c subroutine bcspeval(xget,yget,iselect,fval, > x,nx,y,ny,ilinx,iliny,f,inf3,ier) c real xget,yget real fval(6) real x(nx),y(ny),f(4,4,inf3,ny) c integer iselect(6) integer ilinx,iliny,nx,ny,inf3,ier c c input: c (xget,yget) location where interpolated value is desired c x(1).le.xget.le.x(nx) expected c y(1).le.yget.le.y(ny) expected c c iselect select desired output c c iselect(1)=1 -- want function value (f) itself c iselect(2)=1 -- want df/dx c iselect(3)=1 -- want df/dy c iselect(4)=1 -- want d2f/dx2 c iselect(5)=1 -- want d2f/dy2 c iselect(6)=1 -- want d2f/dxdy c c example: iselect(1)=iselect(2)=iselect(3)=1 c f, df/dx, and df/dy all evaluated c iselect(4)=iselect(5)=iselect(6)=0 c 2nd derivatives not evaluated. c c the number of non zero values iselect(1:6) c determines the number of outputs... c see fval (output) description. 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 iselect(1)=3 -- want 3rd derivatives c iselect(2)=1 for d3f/dx3 c iselect(3)=1 for d3f/dx2dy c iselect(4)=1 for d3f/dxdy2 c iselect(5)=1 for d3f/dy3 c number of non-zero values iselect(2:5) gives no. of outputs c if iselect(1)=4 -- want 4th derivatives c iselect(2)=1 for d4f/dx3dy c iselect(3)=1 for d4f/dx2dy2 c iselect(4)=1 for d4f/dxdy3 c number of non-zero values iselect(2:4) gives no. of outputs c if iselect(1)=5 -- want 5th derivatives c iselect(2)=1 for d5f/dx3dy2 c iselect(3)=1 for d5f/dx2dy3 c number of non-zero values iselect(2:3) gives no. of outputs c if iselect(1)=6 -- want 6th derivatives c d6f/dx3dy3 -- one value is returned. c c x(1...nx) independent coordinate x, strict ascending c y(1...ny) independent coordinate y, strict ascending c c ilinx -- =1: flag that x is linearly spaced c iliny -- =1: flag that y is linearly spaced c c **CAUTION** actual even spacing of x, y is NOT CHECKED HERE! c c c f the function values (at grid points) c and spline coefficients c c evaluation formula: for point x btw x(i) and x(i+1), dx=x-x(i) c and y btw y(j) and y(j+1), dy=y-y(j), c c spline value = c f(1,1,i,j) + dx*f(2,1,i,j) + dx**2*f(3,1,i,j) + dx**3*f(4,1,i,j) c +dy*(f(1,2,i,j) + dx*f(2,2,i,j) + dx**2*f(3,2,i,j) + dx**3*f(4,2,i,j)) c +d2*(f(1,3,i,j) + dx*f(2,3,i,j) + dx**2*f(3,3,i,j) + dx**3*f(4,3,i,j)) c +d3*(f(1,4,i,j) + dx*f(2,4,i,j) + dx**2*f(3,4,i,j) + dx**3*f(4,4,i,j)) c c where d2=dy**2 and d3=dy**3. c c output: c up to 6 elements of fval, ordered as follows: c fval(1)=function value or lowest order derivative requested c fval(2)=next order derivative c etc c the ordering is a subset of the sequence given under the c "iselect" input variable description. c c ier = 0 -- successful completion; = 1 -- an error occurred.

Home | Top |

tcspeval

c tcspeval -- eval tricubic spline function and/or derivatives c subroutine tcspeval(xget,yget,zget,iselect,fval, > x,nx,y,ny,z,nz,ilinx,iliny,ilinz,f,inf4,inf5, > ier) c real xget,yget,zget real fval(10) real x(nx),y(ny),z(nz),f(4,4,4,inf4,inf5,nz) c integer iselect(10) integer ilinx,iliny,ilinz,nx,ny,nz,inf4,inf5,ier c c input: c (xget,yget,zget) location where interpolated value is desired c x(1).le.xget.le.x(nx) expected c y(1).le.yget.le.y(ny) expected c z(1).le.zget.le.z(nz) expected c c iselect select desired output c c iselect(1)=1 -- want function value (f) itself c iselect(2)=1 -- want df/dx c iselect(3)=1 -- want df/dy c iselect(4)=1 -- want df/dz c iselect(5)=1 -- want d2f/dx2 c iselect(6)=1 -- want d2f/dy2 c iselect(7)=1 -- want d2f/dz2 c iselect(8)=1 -- want d2f/dxdy c iselect(9)=1 -- want d2f/dxdz c iselect(10)=1 -- want d2f/dydz c c c example: iselect(1)=iselect(2)=iselect(3)=iselect(4)=1 c f, df/dx, df/dy, and df/dz all evaluated c iselect(5)=iselect(6)=iselect(7)=0 c iselect(8)=iselect(9)=iselect(10)=0 c 2nd derivatives not evaluated. c c see fval (output) description. c c (new dmc Dec 2005 -- higher derivatives available) c iselect(1)=3 --> 3rd derivative, .le.2 diff. in any coordinate c iselect(2:8) select: fxxy fxxz fxyy fxyz fxzz fyyz fyzz c ->note iselect(1)=3, iselect(5)=1 gives fxyz = d3f/dxdydz c iselect(1)=-3 --> 3rd derivative, 3 in one coordinate c iselect(2:4) select: fxxx fyyy fzzz c iselect(1)=4 --> 3rd derivative, .le.2 diff. in any coordinate c iselect(2:7) select: fxxyy fxxyz fxxzz fxyyz fxyzz fyyzz c iselect(1)=-4 --> 3rd derivative, 3 in one coordinate c iselect(2:7) select: fxxxy fxxxz fxyyy fxzzz fyyyz fyzzz c iselect(1)=5 --> 3rd derivative, .le.2 diff. in any coordinate c iselect(2:4) select: fxxyyz fxxyzz fxyyzz c iselect(1)=-5 --> 3rd derivative, 3 in one coordinate c iselect(2:10) select: fxxxyy fxxxyz fxxxzz fxxyyy fxxzzz c fxyyyz fxyzzz fyyyzz fzzzyy c iselect(1)=6 --> 3rd derivative, .le.2 diff. in any coordinate c fxxyyzz c iselect(1)=-6 --> 3rd derivative, 3 in one coordinate c iselect(2:10) select: fxxxyyy fxxxyyz fxxxyzz fxxxyyz c fxxyyyz fxxyzzz fxyyyzz fxyyzzz fyyyzzz c iselect(1)=-7 --> 7th derivative c iselect(2:7) select: fxxxyyyz fxxxyyzz fxxxyzzz c fxxyyyzz fxxyyzzz fxyyyzzz c iselect(1)=-8 --> 8th derivative c iselect(2:4) select: fxxxyyyzz fxxxyyzzz fxxyyyzzz c iselect(1)=-9 --> 9th derivative: fxxxyyyzzz c c x(1...nx) independent coordinate x, strict ascending c y(1...ny) independent coordinate y, strict ascending c z(1...nz) independent coordinate y, strict ascending c c ilinx -- =1: flag that x is linearly spaced c iliny -- =1: flag that y is linearly spaced c ilinz -- =1: flat that z is linearly spaced c c **CAUTION** actual even spacing of x, y, z is NOT CHECKED HERE! c c c f the function values (at grid points) c and spline coefficients c c evaluation formula: for point x btw x(i) and x(i+1), dx=x-x(i) c and y btw y(j) and y(j+1), dy=y-y(j), c and z btw z(k) and z(k+1), dz=z-z(k) c c do m=1,4 c p(m) = c f(1,1,m,i,j,k)+dx*f(2,1,m,i,j,k)+dx**2*f(3,1,m,i,j,k)+dx**3*f(4,1,m,i,j,k) c +dy*( c f(1,2,m,i,j,k)+dx*f(2,2,m,i,j,k)+dx**2*f(3,2,m,i,j,k)+dx**3*f(4,2,m,i,j,k)) c +dy**2*( c f(1,3,m,i,j,k)+dx*f(2,3,m,i,j,k)+dx**2*f(3,3,m,i,j,k)+dx**3*f(4,3,m,i,j,k)) c +dy**3*( c f(1,4,m,i,j,k)+dx*f(2,4,m,i,j,k)+dx**2*f(3,4,m,i,j,k)+dx**3*f(4,4,m,i,j,k)) c enddo c answer = p(1)+dz*p(2)+dz**2*p(3)+dz**3*p(4) c c output: c up to 10 elements of fval, ordered as follows: c fval(1)=function value or lowest order derivative requested c fval(2)=next order derivative c etc c the ordering is a subset of the sequence given under the c "iselect" input variable description; the first M elements c of fval are used, where M = the number of non-zero elements c of iselect. c c ier = 0 -- successful completion; = 1 -- an error occurred. c

Home | Top |

Other_1d_Spline_Routines

Coefficient setup: spline, pspline (Forsythe/Malcolm/Moler style spline coefficients computation) Evaluation: seval, speval (Forsythe/Malcolm/Moler style spline function evaluation)

Home | Top |

spline

C THE CODES (SPLINE & SEVAL) ARE TAKEN FROM: C FORSYTHE,MALCOLM AND MOLER, "COMPUTER METHODS FOR C MATHEMATICAL COMPUTATIONS",PRENTICE-HALL, 1977. C SUBROUTINE SPLINE (N, X, Y, B, C, D) INTEGER N REAL X(N), Y(N), B(N), C(N), D(N) C C THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED C FOR A CUBIC INTERPOLATING SPLINE C C S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3 C C FOR X(I) .LE. X .LE. X(I+1) C C INPUT.. C C N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2) C X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER C Y = THE ORDINATES OF THE KNOTS C C OUTPUT.. C C B, C, D = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE. C

Home | Top |

pspline

C PSPLINE -- modified from SPLINE.FOR -- dmc 3 Oct 1995 C C THE CODES (SPLINE & SEVAL) ARE TAKEN FROM: C FORSYTHE,MALCOLM AND MOLER, "COMPUTER METHODS FOR C MATHEMATICAL COMPUTATIONS",PRENTICE-HALL, 1977. c c PSPLINE -- periodic spline -- adaptation of SPLINE.FOR by D. McCune c October 1995: The function to be splined is taken to have periodic C derivatives (it need not be periodic itself), so that c the interpolating function satisfies c y'(x(1))=y'(x(N)) c y''(x(1))=y''(x(N)) c this results in a fully specified system (no addl BC assumptions c required); it is not a tridiagonal system but a modified NM1xNM1 c tridiagonal system with non-zero offdiagonal corner (1,NM1),(NM1,1) c elements. It remains symmetric & diagonally dominant, and is c solved by Gaussian elimination w/o pivoting. NM1=N-1. c SUBROUTINE PSPLINE (N, X, Y, B, C, D, WK) INTEGER N REAL X(N), Y(N), B(N), C(N), D(N), WK(N) C C THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED C FOR A periodic CUBIC INTERPOLATING SPLINE C C S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3 C C FOR X(I) .LE. X .LE. X(I+1) C C WK(...) is used as a workspace during the Guassian elimination. C it represents the rightmost column of the array C C note B(N),C(N),D(N) give coeffs for a periodic continuation into C the next period, up to X(N)+(X(2)-X(1)) only. C C SEVAL can be used to evaluate the spline, but, it is up to the C SEVAL caller to bring the argument into the normal range C X(1) to X(N). C C INPUT.. C C N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2) C includes a complete period of the data, the last point being C a repeat of the first point. C X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER C X(N)-X(1) is the period of the periodic function represented. C Y = THE ORDINATES OF THE KNOTS C C OUTPUT.. C C B, C, D = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE. C

Home | Top |

seval

(not recommended for MPP use) REAL FUNCTION SEVAL(N, U, X, Y, B, C, D) INTEGER N REAL U, X(N), Y(N), B(N), C(N), D(N) C C THIS SUBROUTINE EVALUATES THE CUBIC SPLINE FUNCTION C C SEVAL = Y(I) + B(I)*(U-X(I)) + C(I)*(U-X(I))**2 + D(I)*(U-X(I))**3 C C WHERE X(I) .LT. U .LT. X(I+1), USING HORNER'S RULE C C IF U .LT. X(1) THEN I = 1 IS USED. C IF U .GE. X(N) THEN I = N IS USED. C C INPUT.. C C N = THE NUMBER OF DATA POINTS C U = THE ABSCISSA AT WHICH THE SPLINE IS TO BE EVALUATED C X,Y = THE ARRAYS OF DATA ABSCISSAS AND ORDINATES C B,C,D = ARRAYS OF SPLINE COEFFICIENTS COMPUTED BY SPLINE C C IF U IS NOT IN THE SAME INTERVAL AS THE PREVIOUS CALL, THEN A C BINARY SEARCH IS PERFORMED TO DETERMINE THE PROPER INTERVAL. C INTEGER I, J, K REAL DX DATA I/1/ IF ( I .GE. N ) I = 1 IF ( U .LT. X(I) ) GO TO 10 IF ( U .LE. X(I+1) ) GO TO 30 C C BINARY SEARCH C 10 I = 1 J = N+1 20 K = (I+J)/2 IF ( U .LT. X(K) ) J = K IF ( U .GE. X(K) ) I = K IF ( J .GT. I+1 ) GO TO 20 C C EVALUATE SPLINE C 30 DX = U - X(I) SEVAL = Y(I) + DX*(B(I) + DX*(C(I) + DX*D(I))) C RETURN END

Home | Top |

speval

(not recommended for MPP use) REAL FUNCTION SPEVAL(N, U, X, Y, B, C, D) INTEGER N REAL U, X(N), Y(N), B(N), C(N), D(N) C C THIS SUBROUTINE EVALUATES THE DERIVATIVE OF THE CUBIC SPLINE C FUNCTION C C C WHERE X(I) .LT. U .LT. X(I+1), USING HORNER'S RULE C C IF U .LT. X(1) THEN I = 1 IS USED. C IF U .GE. X(N) THEN I = N IS USED. C C INPUT.. C C N = THE NUMBER OF DATA POINTS C U = THE ABSCISSA AT WHICH THE SPLINE IS TO BE EVALUATED C X,Y = THE ARRAYS OF DATA ABSCISSAS AND ORDINATES C B,C,D = ARRAYS OF SPLINE COEFFICIENTS COMPUTED BY SPLINE C C IF U IS NOT IN THE SAME INTERVAL AS THE PREVIOUS CALL, THEN A C BINARY SEARCH IS PERFORMED TO DETERMINE THE PROPER INTERVAL. C INTEGER I, J, K REAL DX DATA I/1/ IF ( I .GE. N ) I = 1 IF ( U .LT. X(I) ) GO TO 10 IF ( U .LE. X(I+1) ) GO TO 30 C C BINARY SEARCH C 10 I = 1 J = N+1 20 K = (I+J)/2 IF ( U .LT. X(K) ) J = K IF ( U .GE. X(K) ) I = K IF ( J .GT. I+1 ) GO TO 20 C C EVALUATE SPLINE C 30 DX = U - X(I) SPEVAL = B(I) + DX*(2.*C(I) + 3.*DX*D(I)) RETURN END

Home | Top |

Hermite

Routines for generating Hermite data sets (i.e. derivative evaluation at grid points): *using Akima variation minimization method for derivative evaluation: akherm1, akherm2, akherm3 or Akima with periodic boundary conditions: akherm1p, akherm2p, akherm3p (these routines may be treated as examples of ways to create Hermite data sets) *using calls to user supplied function for derivative evaluation: mkherm1, mkherm2, mkherm3 *using simple numerical differentiation on the given grid: dnherm1, dnherm2, dnherm3 (note: usually better to use the Akima routines) Vectorized Evaluation Routines (recommended): ============================================= given a location vector, do lookup & evaluation: vecherm1 vecherm2 vecherm3 given a grid (different than the original grid), to map the function from the old grid to the new grid: gridherm1 gridherm2 gridherm3 if zone location data is already computed, to do just the vectorized spline evaluation: herm1fcn herm2fcn herm3fcn Scalar Evaluation Routines ========================== herm1ev herm2ev herm3ev

Home | Top |

akherm1p_&_akherm1

subroutine akherm1p(x,nx,fherm,ilinx,ipx,ier) C C Akima subroutine with boundary condition option: C C =>ipx=1 for periodic boundary condition C C =>ipx=2 for user specified boundary condition: C in which case fherm(1,1) and fherm(1,nx) must contain C the derivatives fx(x(1)) and fx(x(nx)), respectively, and these C will not be modified on output. C C =>ipx=0: default boundary conditions C C other arguments as with akherm1. C C input: C integer nx ! array dimensions real x(nx) ! x coordinate array real fherm(0:1,nx) ! data/Hermite array integer ipx ! =1: f' periodic in x C C---------------------------- (or call with default boundary condition...) subroutine akherm1(x,nx,fherm,ilinx,ier) C C create a data set for Hermite interpolation, based on Akima's method C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1] C C 1d routine C C input: C integer nx ! array dimensions real x(nx) ! x coordinate array real fherm(0:1,nx) ! data/Hermite array C C fherm(0,i) = function value f at x(i) **on input** C C fherm(1,i) = derivative df/dx at x(i) **on output** C C addl output: C ilinx=1 if x is "evenly spaced" ier=0 if no errors C C ** x must be strict ascending ** C C ier =0 on out if there is no error. C

Home | Top |

akherm2p_&_akherm2

subroutine akherm2p(x,nx,y,ny,fherm,nf2,ilinx,iliny,ipx,ipy,ier) C C create a data set for Hermite interpolation, based on Akima's method C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1] C C with independently settable boundary condition options: C ipx or ipy =0 -- default, boundary conditions from divided diffs C ipx or ipy =1 -- periodic boundary condition C ipx or ipy =2 -- user supplied df/dx or df/dy C input: C integer nx,ny,nf1 ! array dimensions real x(nx) ! x coordinate array real y(ny) ! y coordinate array real fherm(0:3,nf2,ny) ! data/Hermite array C integer ipx ! =1 if df/dx periodic in x integer ipy ! =1 if df/dy periodic in y C C fherm(0,1:nx,1:ny) supplied by user; this routine computes the C rest of the elements, but note: C C if ipx=2: fherm(1,1,1:ny) & fherm(1,nx,1:ny) = INPUT df/dx BCs C if ipy=2: fherm(2,1:nx,1) & fherm(2,1:nx,ny) = INPUT df/dy BCs C C on output, at all grid points (i,j) covering (1:nx,1:ny): C fherm1(1,i,j) -- df/dx at the grid point C fherm1(2,i,j) -- df/dy at the grid point C fherm1(3,i,j) -- d2f/dxdy at the grid point C C--------------------- (or call with default boundary condition...) subroutine akherm2(x,nx,y,ny,fherm,nf2,ilinx,iliny,ier) C C create a data set for Hermite interpolation, based on Akima's method C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1] C C input: C integer nx,ny,nf1 ! array dimensions real x(nx) ! x coordinate array real y(ny) ! y coordinate array real fherm(0:3,nf2,ny) ! data/Hermite array C C fherm(0,i,j) = function value f at x(i),y(j) **on input** C C fherm(1,i,j) = derivative df/dx at x(i),y(j) **on output** C fherm(2,i,j) = derivative df/dy at x(i),y(j) **on output** C fherm(3,i,j) = derivative d2f/dxdy at x(i),y(j) **on output** C C addl output: C ilinx=1 if x axis is evenly spaced C iliny=1 if y axis is evenly spaced C ier=0 if no error: C x, y must both be strict ascending C nf2.ge.nx is required. C

Home | Top |

akherm3p_&akherm3

C---------------------------- subroutine akherm3p(x,nx,y,ny,z,nz,fherm,nf2,nf3, > ilinx,iliny,ilinz,ipx,ipy,ipz,ier) C C create a data set for Hermite interpolation, based on Akima's method C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1] C C 3d routine -- C with independently settable boundary condition options: C ipx|ipy|ipz =0 -- default, boundary conditions from divided diffs C ipx|ipy|ipz =1 -- periodic boundary condition C ipx|ipy|ipz =2 -- user supplied df/dx|df/dy|df/dz C C input: C C nf2.get.nx and nf3.ge.ny expected! C integer nx,ny,nz,nf2,nf3 ! array dimensions real x(nx) ! x coordinate array real y(ny) ! y coordinate array real z(nz) ! z coordinate array real fherm(0:7,nf2,nf3,nz) ! data/Hermite array C integer ipx ! =0: dflt; =1: periodic df/dx integer ipy ! =0: dflt; =1: periodic df/dy integer ipz ! =0: dflt; =1: periodic df/dz C C if ipx=2: fherm(1,1,1:ny,1:nz) & fherm(1,nx,1:ny,1:nz) = INPUT df/dx BCs C if ipy=2: fherm(2,1:nx,1,1:nz) & fherm(2,1:nx,ny,1:nz) = INPUT df/dy BCs C if ipz=2: fherm(3,1:nx,1:ny,1) & fherm(3,1:nx,1:ny,nz) = INPUT df/dz BCs C C for the rest of the grid points... C C fherm(0,i,j,k) = function value f at x(i),y(j),z(k) **on input** C C fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k) **on output** C fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k) **on output** C fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k) **on output** C fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k) **on output** C fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k) **on output** C fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k) **on output** C fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k) **on output** C C additional outputs: C C ilinx = 1 if x is evenly spaced (approximately) C iliny = 1 if y is evenly spaced (approximately) C ilinz = 1 if z is evenly spaced (approximately) C C--------------------- (or call with default boundary condition...) subroutine akherm3(x,nx,y,ny,z,nz,fherm,nf2,nf3, > ilinx,iliny,ilinz,ier) C C create a data set for Hermite interpolation, based on Akima's method C [Hiroshi Akima, Communications of the ACM, Jan 1974, Vol. 17 No. 1] C C 3d routine C C input: C C nf2.get.nx and nf3.ge.ny expected! C integer nx,ny,nz,nf2,nf3 ! array dimensions real x(nx) ! x coordinate array real y(ny) ! y coordinate array real z(nz) ! z coordinate array real fherm(0:7,nf2,nf3,nz) ! data/Hermite array C C fherm(0,i,j,k) = function value f at x(i),y(j),z(k) **on input** C C fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k) **on output** C fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k) **on output** C fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k) **on output** C fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k) **on output** C fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k) **on output** C fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k) **on output** C fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k) **on output** C C additional outputs: C C ilinx = 1 if x is evenly spaced (approximately) C iliny = 1 if y is evenly spaced (approximately) C ilinz = 1 if z is evenly spaced (approximately) C C ier = 0 if there are no errors C C note possible errors: x y or z NOT strict ascending C nf2.lt.nx .or. nf3.lt.ny C

Home | Top |

mkherm1

subroutine mkherm1(fun,x,nx,fherm) C C create a data set for Hermite interpolation, from evaluation of C function and derivatives. function of 1 indep. coordinate. C C input: C external fun ! passed real function(x,dfdx) real x(nx) ! x coordinate array C C the passed function fun must have the interface: C C real function <name>(x,dfdx) C where x is input, the function returns the function value and the C argument dfdx returns as output the function derivative at x. C C output: C real fherm(0:1,nx) ! function data & derivative C C fherm(0,j) = function value f at x(j) C fherm(1,j) = derivative df/dx at x(j) C C---------------------------- C do ix=1,nx fherm(0,ix)=fun(x(ix),dfdx) fherm(1,ix)=dfdx enddo C return end

Home | Top |

mkherm2

subroutine mkherm2(fun,x,nx,y,ny,fherm) C C create a data set for Hermite interpolation, from evaluation of C function and derivatives. function of 2 indep. coordinates. C C input: C external fun ! passed real function(x,y) real x(nx) ! x coordinate array real y(ny) ! y coordinate array C C the passed function fun must have the interface: C C real function <name>(x,y,dfdx,dfdy,d2fdxdy) C where x,y are input, the function returns the function value, C and the arguments dfdx and dfdy return as output the function C derivative at the point (x,y). C C output: C real fherm(0:3,nx,ny) ! function data & derivatives C C fherm(0,i,j) = function value f at x(i),y(j) C fherm(1,i,j) = derivative df/dx at x(i),y(j) C fherm(2,i,j) = derivative df/dy at x(i),y(j) C fherm(3,i,j) = derivative d2f/dxdy at x(i),y(j) C C---------------------------- C do ix=1,nx do iy=1,ny fherm(0,ix,iy)=fun(x(ix),y(iy),dfdx,dfdy,d2fdxdy) fherm(1,ix,iy)=dfdx fherm(2,ix,iy)=dfdy fherm(3,ix,iy)=d2fdxdy enddo enddo C return end

Home | Top |

mkherm3

subroutine mkherm3(fun,x,nx,y,ny,z,nz,fherm) C C create a data set for Hermite interpolation, from evaluation of C function and derivatives. Function of 3 indep. coordinates. C C input: C external fun ! passed real function(x,y,z) real x(nx) ! x coordinate array (1st dim) real y(ny) ! y coordinate array (2nd dim) real z(nz) ! z coordinate array (3rd dim) C C the passed function fun must have the interface: C C real function <name>(x,y,z,dfdx,dfdy,dfdz, C d2fdxdy,d2fdxdz,d2fdydz,d3fdxdydz) C where x,y,z are input, the function returns the function value, C and the arguments dfdx, dfdy and dfdz return as output the function C derivative at the point (x,y,z). C C output: C real fherm(0:7,nx,ny,nz) ! function data & derivatives C C fherm(0,i,j,k) = function value f at x(i),y(j),z(k) C fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k) C fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k) C fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k) C fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k) C fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k) C fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k) C fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k) C C---------------------------- C do iz=1,nz do iy=1,ny do ix=1,nx fherm(0,ix,iy,iz)= > fun(x(ix),y(iy),z(iz),dfdx,dfdy,dfdz, > d2fdxdy,d2fdxdz,d2fdydz,d3fdxyz) fherm(1,ix,iy,iz)=dfdx fherm(2,ix,iy,iz)=dfdy fherm(3,ix,iy,iz)=dfdz fherm(4,ix,iy,iz)=d2fdxdy fherm(5,ix,iy,iz)=d2fdxdz fherm(6,ix,iy,iz)=d2fdydz fherm(7,ix,iy,iz)=d3fdxyz enddo enddo enddo C return end

Home | Top |

dnherm1

subroutine dnherm1(x,nx,fherm,ilinx,ier) C C create a data set for Hermite interpolation, based on simple C numerical differentiation using the given grid. C C 1d routine C C input: C integer nx ! array dimensions real x(nx) ! x coordinate array real fherm(0:1,nx) ! data/Hermite array C C fherm(0,i) = function value f at x(i) **on input** C C fherm(1,i) = derivative df/dx at x(i) **on output** C C addl output: C ilinx=1 if x is "evenly spaced" ier=0 if no errors C C ** x must be strict ascending ** C C---------------------------- C ztol=1.0e-3 call splinck(x,nx,ilinx,ztol,ier) if(ier.ne.0) then write(6,*) '?dnherm1: x axis not strict ascending.' return endif C do ix=1,nx c c x div. diffs in vicinity c ixp=min(nx,ix+1) ixm=max(1,ix-1) zd=(fherm(0,ixp)-fherm(0,ixm))/(x(ixp)-x(ixm)) c fherm(1,ix)=zd enddo C return end

Home | Top |

dnherm2

subroutine dnherm2(x,nx,y,ny,fherm,nf2,ilinx,iliny,ier) C C create a data set for Hermite interpolation, based on simple C numerical differentiation using the given grid. The grid should C be "approximately" evenly spaced. C C 2d routine C C input: C C nf2.gt.nx expected! C integer nx,ny,nf2 ! array dimensions real x(nx) ! x coordinate array real y(ny) ! y coordinate array real fherm(0:3,nf2,ny) ! data/Hermite array C C fherm(0,i,j) = function value f at x(i),y(j) **on input** C C fherm(1,i,j) = derivative df/dx at x(i),y(j) **on output** C fherm(2,i,j) = derivative df/dy at x(i),y(j) **on output** C fherm(3,i,j) = derivative d2f/dxdy at x(i),y(j) **on output** C C addl output: C ilinx=1 if x axis is evenly spaced C iliny=1 if y axis is evenly spaced C ier=0 if no error: C x, y must both be strict ascending C nf2.ge.nx is required. C C---------------------------- c ier=0 c call splinck(x,nx,ilinx,1.0e-3,ierx) if(ierx.ne.0) ier=2 c if(ier.eq.2) then write(6,'('' ?dnherm2: x axis not strict ascending'')') endif c call splinck(y,ny,iliny,1.0e-3,iery) if(iery.ne.0) ier=3 c if(ier.eq.3) then write(6,'('' ?dnherm2: y axis not strict ascending'')') endif c if(nf2.lt.nx) then ier=4 write(6,*) '?dnherm2: fherm array dimension too small.' endif C if(ier.ne.0) return C do iy=1,ny c iyp=min(ny,iy+1) iym=max(1,iy-1) c do ix=1,nx c ixp=min(nx,ix+1) ixm=max(1,ix-1) c c x diffs in vicinity c zd=(fherm(0,ixp,iy)-fherm(0,ixm,iy))/(x(ixp)-x(ixm)) c fherm(1,ix,iy)=zd c c y diffs in vicinity c zd=(fherm(0,ix,iyp)-fherm(0,ix,iym))/(y(iyp)-y(iym)) c fherm(2,ix,iy)=zd c c xy cross derivs (except at corners, iedge=2) c fherm(3,ix,iy)=(fherm(0,ixp,iyp)-fherm(0,ixm,iyp) > -fherm(0,ixp,iym)+fherm(0,ixm,iym))/ > ((x(ixp)-x(ixm))*(y(iyp)-y(iym))) c enddo enddo C return end

Home | Top |

dnherm3

subroutine dnherm3(x,nx,y,ny,z,nz,fherm,nf2,nf3, > ilinx,iliny,ilinz,ier) C C create a data set for Hermite interpolation, based on simple C numerical differentiation using the given grid. The grid should C be "approximately" evenly spaced. C C 3d routine C C input: C C nf2.get.nx and nf3.ge.ny expected! C integer nx,ny,nz,nf2,nf3 ! array dimensions real x(nx) ! x coordinate array real y(ny) ! y coordinate array real z(nz) ! z coordinate array real fherm(0:7,nf2,nf3,nz) ! data/Hermite array C C fherm(0,i,j,k) = function value f at x(i),y(j),z(k) **on input** C C fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k) **on output** C fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k) **on output** C fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k) **on output** C fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k) **on output** C fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k) **on output** C fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k) **on output** C fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k) **on output** C C additional outputs: C C ilinx = 1 if x is evenly spaced (approximately) C iliny = 1 if y is evenly spaced (approximately) C ilinz = 1 if z is evenly spaced (approximately) C C ier = 0 if there are no errors C C note possible errors: x y or z NOT strict ascending C nf2.lt.nx .or. nf3.lt.ny C C---------------------------- C C C error check C ier=0 c call splinck(x,nx,ilinx,1.0e-3,ierx) if(ierx.ne.0) ier=2 c if(ier.eq.2) then write(6,'('' ?dnherm3: x axis not strict ascending'')') endif c call splinck(y,ny,iliny,1.0e-3,iery) if(iery.ne.0) ier=3 c if(ier.eq.3) then write(6,'('' ?dnherm3: y axis not strict ascending'')') endif c call splinck(z,nz,ilinz,1.0e-3,ierz) if(ierz.ne.0) ier=4 c if(ier.eq.4) then write(6,'('' ?dnherm3: z axis not strict ascending'')') endif c if(nf2.lt.nx) then ier=5 write(6,*) '?dnherm3: fherm (x) array dimension too small.' endif c if(nf3.lt.ny) then ier=6 write(6,*) '?dnherm3: fherm (y) array dimension too small.' endif C if(ier.ne.0) return C C error check OK C do iz=1,nz c izp=min(nz,iz+1) izm=max(1,iz-1) c do iy=1,ny c iyp=min(ny,iy+1) iym=max(1,iy-1) c do ix=1,nx c ixp=min(nx,ix+1) ixm=max(1,ix-1) c c x diffs in vicinity c zd=(fherm(0,ixp,iy,iz)-fherm(0,ixm,iy,iz))/ > (x(ixp)-x(ixm)) c fherm(1,ix,iy,iz)=zd c c y diffs in vicinity c zd=(fherm(0,ix,iyp,iz)-fherm(0,ix,iym,iz))/ > (y(iyp)-y(iym)) c fherm(2,ix,iy,iz)=zd c c z diffs in vicinity c zd=(fherm(0,ix,iy,izp)-fherm(0,ix,iy,izm))/ > (z(izp)-z(izm)) c fherm(3,ix,iy,iz)=zd c c xy cross derivs c fherm(4,ix,iy,iz)= > (fherm(0,ixp,iyp,iz)-fherm(0,ixm,iyp,iz) > -fherm(0,ixp,iym,iz)+fherm(0,ixm,iym,iz))/ > ((x(ixp)-x(ixm))*(y(iyp)-y(iym))) c c xz cross derivs c fherm(5,ix,iy,iz)= > (fherm(0,ixp,iy,izp)-fherm(0,ixm,iy,izp) > -fherm(0,ixp,iy,izm)+fherm(0,ixm,iy,izm))/ > ((x(ixp)-x(ixm))*(z(izp)-z(izm))) c c yz cross derivs c fherm(6,ix,iy,iz)= > (fherm(0,ix,iyp,izp)-fherm(0,ix,iym,izp) > -fherm(0,ix,iyp,izm)+fherm(0,ix,iym,izm))/ > ((y(iyp)-y(iym))*(z(izp)-z(izm))) c c xyz cross deriv c fherm(7,ix,iy,iz)= > ((fherm(0,ixp,iyp,izp)-fherm(0,ixp,iym,izp) > -fherm(0,ixp,iyp,izm)+fherm(0,ixp,iym,izm))- > (fherm(0,ixm,iyp,izp)-fherm(0,ixm,iym,izp) > -fherm(0,ixm,iyp,izm)+fherm(0,ixm,iym,izm)))/ > ((x(ixp)-x(ixm))*(y(iyp)-y(iym))*(z(izp)-z(izm))) c enddo enddo enddo C return end

Home | Top |

vecherm1

subroutine vecherm1(ict,ivec,xvec,ivd,fval,nx,xpkg,fspl, > iwarn,ier) c c vectorized hermite evaluation routine -- 1d c 1. call vectorized zone lookup routine c 2. call vectorized hermite 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 integer ivec ! vector dimensioning c c ivec-- number of vector pts (hermite 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 hermite x grid real xpkg(nx,4) ! x grid "package" (cf genxpkg) real fspl(2,nx) ! Hermite coefficients (cf herm1ev) 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

Home | Top |

vecherm2

subroutine vecherm2(ict,ivec,xvec,yvec,ivd,fval, > nx,xpkg,ny,ypkg,fspl,inf2, > iwarn,ier) c c vectorized hermite evaluation routine -- 2d c 1. call vectorized zone lookup routine c 2. call vectorized hermite 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/dxdy (don't evaluate if ict(4)=0) c integer ivec ! vector dimensioning c c ivec-- number of vector pts (hermite 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 hermite 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) ! Hermite coefficients, cf herm2ev 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

Home | Top |

vecherm3

subroutine vecherm3(ict,ivec,xvec,yvec,zvec,ivd,fval, > nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf2,inf3, > iwarn,ier) c c vectorized hermite evaluation routine -- c 1. call vectorized zone lookup routine c 2. call vectorized hermite evaluation routine c c-------------------------- c input: integer ict(8) ! 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/dxdy (don't evaluate if ict(5)=0) c ict(6)=1 for d2f/dxdz (don't evaluate if ict(6)=0) c ict(7)=1 for d2f/dydz (don't evaluate if ict(7)=0) c ict(8)=1 for d3f/dxdydz (don't evaluate if ict(8)=0) c integer ivec ! vector dimensioning c c ivec-- number of vector pts (hermite 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 hermite 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) ! Hermite coefficients cf herm3ev 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

Home | Top |

gridherm1

subroutine gridherm1(x_newgrid,nx_new,f_new,nx,xpkg,fspl, > iwarn,ier) c c regrid a hermite 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(2) c data ict/1,0/ c c-------------------------------------------- c call vecherm1(ict,nx_new,x_newgrid,nx_new,f_new,nx,xpkg,fspl, > iwarn,ier) c return end

Home | Top |

gridherm2

subroutine gridherm2( > x_newgrid,nx_new,y_newgrid,ny_new,f_new,if1, > nx,xpkg,ny,ypkg,fspl,inf2,iwarn,ier) c c regrid a hermite 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 inf2 ! array dimension real fspl(0:4,inf2,ny) ! hermite 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

Home | Top |

gridherm3

subroutine gridherm3( > x_newgrid,nx_new,y_newgrid,ny_new,z_newgrid,nz_new, > f_new,if1,if2, > nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf2,inf3,iwarn,ier) c c regrid a hermite 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 inf2,inf3 ! array dimensions real fspl(0:7,inf2,inf3,nz) ! hermite 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

Home | Top |

herm1fcn

C evaluate C1 cubic Hermite function interpolation -- 1d fcn C --vectorized-- dmc 10 Feb 1999 C subroutine herm1fcn(ict,ivec,ivecd, > fval,ii,xparam,hx,hxi,fin,nx) C C input: C integer ict(2) ! requested output control integer ivec ! vector length integer ivecd ! vector dimension (1st dim of fval) C integer ii(ivec) ! target cells real xparam(ivec) ! normalized displacements in cells 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) ! Hermite dataset C C output: C real fval(ivecd,2) ! interpolation results C C for detailed description of fin, ict and fval see subroutine C herm1ev 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 C output array fval has a vector ** 1st dimension ** whose C size must be given as a separate argument; ivecd.ge.ivec C expected! C C to use this routine in scalar mode, pass in ivec=ivecd=1 C

Home | Top |

herm2fcn

C evaluate C1 cubic Hermite function interpolation -- 2d fcn C --vectorized-- dmc 10 Feb 1999 C subroutine herm2fcn(ict,ivec,ivecd, > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi, > fin,inf2,ny) C integer ict(4) ! requested output selector 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 "herm2ev") C ! f,fx,fy,fxy at each grid point. C real fval(ivecd,4) ! output returned C ! function value, derivatives... C C for detailed description of fin, ict and fval see subroutine C herm2ev 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 C

Home | Top |

herm3fcn

C evaluate C1 cubic Hermite function interpolation -- 3d fcn C --vectorized-- dmc 10 Feb 1999 C subroutine herm3fcn(ict,ivec,ivecd, > fval,ii,jj,kk,xparam,yparam,zparam, > hx,hxi,hy,hyi,hz,hzi, > fin,inf2,inf3,nz) C integer ict(8) ! requested output selector 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 "herm3ev") C ! f,fx,fy,fz,fxy,fxz,fyz,fxyz at each C ! grid point C real fval(ivecd,8) ! output returned C ! function value, derivatives... C C for detailed description of fin, ict and fval see subroutine herm3ev C 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

Home | Top |

herm1ev

subroutine herm1ev(xget,x,nx,ilinx,f,ict,fval,ier) C C evaluate a 1d cubic Hermite interpolant -- this is C1. C C input arguments: C ================ C 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 contents of f: C C f(0,i) = f @ x(i) C f(1,i) = df/dx @ x(i) C integer ict(2) ! 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 C output arguments: C ================= C real fval(2) ! 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 C examples: C on input ict = [1,1] C on output fval= [f,df/dx] C C on input ict = [1,0] C on output fval= [f] ... element 2 never referenced C C on input ict = [0,1] C on output fval= [df/dx] ... element 2 never referenced C C ier -- completion code: 0 means OK

Home | Top |

herm2ev

subroutine herm2ev(xget,yget,x,nx,y,ny,ilinx,iliny, > f,inf2,ict,fval,ier) C C evaluate a 2d cubic Hermite interpolant on a rectilinear C grid -- this is C1 in both directions. C C this subroutine calls two subroutines: C herm2xy -- find cell containing (xget,yget) C herm2fcn -- evaluate interpolant function and (optionally) C the derivatives C C input arguments: C ================ C 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) = df/dx @ x(i),y(j) C f(2,i,j) = df/dy @ x(i),y(j) C f(3,i,j) = d2f/dxdy @ x(i),y(j) 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/dxdy (0, don't) C C output arguments: C ================= C real fval(4) ! 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 C examples: C on input ict = [1,1,1,1] C on output fval= [f,df/dx,df/dy,d2f/dxdy] C C on input ict = [1,0,0,0] C on output fval= [f] ... elements 2 & 3 & 4 never referenced C C on input ict = [0,1,1,0] C on output fval= [df/dx,df/dy] ... element 3 & 4 never referenced C C on input ict = [0,0,1,0] C on output fval= [df/dy] ... elements 2 & 3 & 4 never referenced. C C ier -- completion code: 0 means OK

Home | Top |

herm3ev

subroutine herm3ev(xget,yget,zget,x,nx,y,ny,z,nz, > ilinx,iliny,ilinz, > f,inf2,inf3,ict,fval,ier) C C evaluate a 3d cubic Hermite interpolant on a rectilinear C grid -- this is C1 in all directions. C C this subroutine calls two subroutines: C herm3xyz -- find cell containing (xget,yget,zget) C herm3fcn -- evaluate interpolant function and (optionally) C the derivatives 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) = df/dx @ x(i),y(j),z(k) C f(2,i,j,k) = df/dy @ x(i),y(j),z(k) C f(3,i,j,k) = df/dz @ x(i),y(j),z(k) C f(4,i,j,k) = d2f/dxdy @ x(i),y(j),z(k) C f(5,i,j,k) = d2f/dxdz @ x(i),y(j),z(k) C f(6,i,j,k) = d2f/dydz @ x(i),y(j),z(k) C f(7,i,j,k) = d3f/dxdydz @ x(i),y(j),z(k) C integer ict(8) ! 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/dxdy (0, don't) C ict(6)=1 -- return d2f/dxdz (0, don't) C ict(7)=1 -- return d2f/dydz (0, don't) C ict(8)=1 -- return d3f/dxdydz (0, don't) C C output arguments: C ================= C real fval(8) ! 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-8) receive 5th thru 8th outputs (if req. by ict(...) spec) C C examples: C on input ict = [1,1,1,1,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] C on output fval= [f] ... elements 2-8 never referenced C C on input ict = [0,1,1,0,0,0,0,0] C on output fval= [df/dx,df/dy] ... elements 3-8 never referenced C C on input ict = [0,0,0,0,1,0,0,0] C on output fval= [d2f/dxdy] ... elements 2-8 never referenced. C C ier -- completion code: 0 means OK

Home | Top |

Piecewise_Linear

The pspline library also contains routines for 1d, 2d, and 3d piecewise linear interpolation. Piecewise linear interpolation is fast (most of the work is in zone lookup), and the user need only supply data at the grid points; no spline or Hermite coefficients are required. But, the interpolant, while continuous, is not differentiable across grid zone boundaries. There are no setup routines for the data itself (use of genxpkg is still recommended for the grids). Vectorized Evaluation Routines (recommended): ============================================= given a location vector, do lookup & evaluation: vecpc1 vecpc2 vecpc3 if zone location data is already computed, to do just the vectorized interpolation evaluation: pc1fcn pc2fcn pc3fcn Scalar Evaluation Routines ========================== pc1ev pc2ev pc3ev

Home | Top |

vecpc1

subroutine vecpc1(ict,ivec,xvec,ivd,fval,nx,xpkg,fspl, > iwarn,ier) c c vectorized piecewise linear evaluation routine -- 1d c 1. call vectorized zone lookup routine c 2. call vectorized piecewise linear evaluation routine c c note -- derivatives are *not* continuous 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 integer ivec ! vector dimensioning c c ivec-- number of vector pts (piecewise linear 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 x grid real xpkg(nx,4) ! x grid "package" (cf genxpkg) real fspl(nx) ! Piecewise Linear function 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

Home | Top |

vecpc2

subroutine vecpc2(ict,ivec,xvec,yvec,ivd,fval, > nx,xpkg,ny,ypkg,fspl,inf1, > iwarn,ier) c c vectorized piecewise linear evaluation routine -- 2d c 1. call vectorized zone lookup routine c 2. call vectorized piecewise linear 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/dxdy (don't evaluate if ict(4)=0) c c note -- derivatives are *not* continuous! c integer ivec ! vector dimensioning c c ivec-- number of vector pts (piecewise linear 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 axis grids real xpkg(nx,4) ! x grid "package" (cf genxpkg) real ypkg(ny,4) ! y grid "package" (cf genxpkg) integer inf1 ! fspl 3rd array dimension, .ge.nx real fspl(inf1,ny) ! Piecewise Linear function 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

Home | Top |

vecpc3

subroutine vecpc3(ict,ivec,xvec,yvec,zvec,ivd,fval, > nx,xpkg,ny,ypkg,nz,zpkg,fspl,inf1,inf2, > iwarn,ier) c c vectorized piecewise linear evaluation routine -- c 1. call vectorized zone lookup routine c 2. call vectorized piecewise linear evaluation routine c c-------------------------- c input: integer ict(8) ! 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/dxdy (don't evaluate if ict(5)=0) c ict(6)=1 for d2f/dxdz (don't evaluate if ict(6)=0) c ict(7)=1 for d2f/dydz (don't evaluate if ict(7)=0) c ict(8)=1 for d3f/dxdydz (don't evaluate if ict(8)=0) c c note -- derivatives are *not* continuous! c integer ivec ! vector dimensioning c c ivec-- number of vector pts (piecewise linear 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 axis 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 inf1 ! fspl 1st array dimension, .ge.nx integer inf2 ! fspl 2nd array dimension, .ge.ny real fspl(0:7,inf1,inf2,nz) ! Piecewise Linear function array 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

Home | Top |

pc1fcn

C--------------------------------------------------------------------- C evaluate C0 piecewise linear function interpolation -- 1d fcn C --vectorized-- dmc 10 Feb 1999 C subroutine pc1fcn(ict,ivec,ivecd, > fval,ii,xparam,hx,hxi,fin,nx) C C input: C integer ict(2) ! requested output control integer ivec ! vector length integer ivecd ! vector dimension (1st dim of fval) C integer ii(ivec) ! target cells real xparam(ivec) ! normalized displacements in cells C real hx(ivec) ! grid spacing, and real hxi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i)) C real fin(nx) ! the data C C output: C real fval(ivecd,2) ! interpolation results C C for detailed description of fin, ict and fval see subroutine C pc1ev 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 C output array fval has a vector ** 1st dimension ** whose C size must be given as a separate argument; ivecd.ge.ivec C expected! C C to use this routine in scalar mode, pass in ivec=ivecd=1 C

Home | Top |

pc2fcn

C--------------------------------------------------------------------- C evaluate piecewise bilinear function interpolation -- 2d fcn C --vectorized-- dmc 10 Feb 1999 C subroutine pc2fcn(ict,ivec,ivecd, > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi, > fin,inf2,ny) C integer ict(4) ! 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(inf2,ny) ! interpolant data (cf "pc2ev") C real fval(ivecd,4) ! output returned C C for detailed description of fin, ict and fval see subroutine C pc2ev 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

Home | Top |

pc3fcn

C--------------------------------------------------------------------- C evaluate trilinear function interpolation -- 3d fcn C --vectorized-- dmc 10 Feb 1999 C subroutine pc3fcn(ict,ivec,ivecd, > fval,ii,jj,kk,xparam,yparam,zparam, > hx,hxi,hy,hyi,hz,hzi, > fin,inf2,inf3,nz) C integer ict(8) ! 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(inf2,inf3,nz) ! interpolant data (cf "pc3ev") C real fval(ivecd,8) ! output returned C C for detailed description of fin, ict and fval see subroutine pc3ev C 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 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

Home | Top |

pc1ev

subroutine pc1ev(xget,x,nx,ilinx,f,ict,fval,ier) C C evaluate a 1d piecewise linear interpolant -- this is C0. C ...a derivative can be evaluated but it is not continuous. C C this subroutine calls two subroutines: C herm1x -- find cell containing (xget,yget) C pc1fcn -- evaluate interpolant function and (optionally) derivatives C C input arguments: C ================ C real xget ! target of this interpolation real x(nx) ! ordered x grid integer ilinx ! ilinx=1 => assume x evenly spaced C real f(nx) ! function data C C contents of f: C C f(i) = f @ x(i) C integer ict(2) ! 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 C output arguments: C ================= C real fval(2) ! 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 C examples: C on input ict = [1,1] C on output fval= [f,df/dx] C C on input ict = [1,0] C on output fval= [f] ... element 2 never referenced C C on input ict = [0,1] C on output fval= [df/dx] ... element 2 never referenced C C ier -- completion code: 0 means OK

Home | Top |

pc2ev

subroutine pc2ev(xget,yget,x,nx,y,ny,ilinx,iliny, > f,inf2,ict,fval,ier) C C evaluate a piecewise bilinear interpolant on a rectilinear C grid -- this is C1 in both directions. C C this subroutine calls two subroutines: C herm2xy -- find cell containing (xget,yget) C pc2fcn -- evaluate interpolant function and (optionally) derivatives C C input arguments: C ================ C 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(inf2,ny) ! function data C C f 2nd dimension inf2 must be .ge. nx C contents of f: C C f(i,j) = f @ x(i),y(j) 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/dxdy (0, don't) C C output arguments: C ================= C real fval(4) ! 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 C examples: C on input ict = [1,1,1,1] C on output fval= [f,df/dx,df/dy,d2f/dxdy] C C on input ict = [1,0,0,0] C on output fval= [f] ... elements 2 & 3 & 4 never referenced C C on input ict = [0,1,1,0] C on output fval= [df/dx,df/dy] ... element 3 & 4 never referenced C C on input ict = [0,0,1,0] C on output fval= [df/dy] ... elements 2 & 3 & 4 never referenced. C C ier -- completion code: 0 means OK

Home | Top |

pc3ev

subroutine pc3ev(xget,yget,zget,x,nx,y,ny,z,nz, > ilinx,iliny,ilinz, > f,inf2,inf3,ict,fval,ier) C C evaluate a trilinear interpolant on a rectilinear grid C derivatives are available, but, not continuous across grid planes. C C this subroutine calls two subroutines: C herm3xyz -- find cell containing (xget,yget,zget) C pc3fcn -- evaluate interpolant function and (optionally) derivatives 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(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(i,j,k) = f @ x(i),y(j),z(k) C integer ict(8) ! 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/dxdy (0, don't) C ict(6)=1 -- return d2f/dxdz (0, don't) C ict(7)=1 -- return d2f/dydz (0, don't) C ict(8)=1 -- return d3f/dxdydz (0, don't) C C output arguments: C ================= C real fval(8) ! 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-8) receive 5th thru 8th outputs (if required by ict(...) spec) C C examples: C on input ict = [1,1,1,1,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] C on output fval= [f] ... elements 2-8 never referenced C C on input ict = [0,1,1,0,0,0,0,0] C on output fval= [df/dx,df/dy] ... elements 3-8 never referenced C C on input ict = [0,0,0,0,1,0,0,0] C on output fval= [d2f/dxdy] ... elements 2-8 never referenced. C C ier -- completion code: 0 means OK

Home | Top |

v_spline

The v_spline subroutine, by Wayne Houlberg, is used throughout the pspline library to actually solve for 1d spline coefficients. Caution -- the convention for representation of non-compact spline coefficients is slightly different between v_spline and the rest of pspline. In v_spline, the convention is f(i,j) = (i-1)'th derivative of the spline at the left hand side of interval j and the formulae for evaluation, for x(j).le.x.le.x(j+1), is delx=x-x(j) s1(x) = f(1,j)+ delx*(f(2,j)+delx*(1.0/2.0)*(f(3,j)+delx*(1.0/3.0)*f(4,j))) s1'(x) = f(2,j)+delx*(f(3,j)+delx*(1.0/2.0)*f(4,j)) s1''(x) = f(3,j)+delx*f(4,j) s1'''(x) = f(4,j) But, in the explicit spline coefficients used by pspline (e.g. the 1d routine cspline), the f(i,j) is the coefficient of the (i-1)'th power term in the cubic polynomial, i.e. the formula is s2(x) = f(1,j)+ delx*(f(2,j)+delx*(f(3,j)+delx*f(4,j))) s2'(x) = f(2,j)+2.0*delx*(f(3,j)+1.5*delx*f(4,j)) s2''(x) = 2.0*f(3,j) + 6.0*delx*f(4,j) s2'''(x) = 6.0*f(4,j) i.e. pspline's f(3,j) = (1.0/2.0) * v_spline's f(3,j) pspline's f(4,j) = (1.0/6.0) * v_spline's f(4,j) Bearing this in mind, users who want access to v_spline directly may do so. The interface is as follows: SUBROUTINE V_SPLINE(k_bc1,k_bcn,n,x,f,wk) !*********************************************************************** !V_SPLINE evaluates the coefficients for a 1d cubic interpolating spline !References: ! Forsythe, Malcolm, Moler, Computer Methods for Mathematical ! Computations, Prentice-Hall, 1977, p.76 ! Engeln-Muellges, Uhlig, Numerical Algorithms with Fortran, Springer, ! 1996, p.251 ! W.A.Houlberg, D.McCune 3/2000 !Input: ! k_bc1-option for BC at x(1) ! =-1 periodic, ignore k_bcn ! =0 not-a-knot ! =1 s'(x1) = input value of f(2,1) ! =2 s''(x1) = input value of f(3,1) ! =3 s'(x1) = 0.0 ! =4 s''(x1) = 0.0 ! =5 match first derivative to first 2 points ! =6 match second derivative to first 3 points ! =7 match third derivative to first 4 points ! =else use not-a-knot ! k_bcn-option for boundary condition at x(n) ! =0 not-a-knot ! =1 s'(x1) = input value of f(2,1) ! =2 s''(x1) = input value of f(3,1) ! =3 s'(x1) = 0.0 ! =4 s''(x1) = 0.0 ! =5 match first derivative to first 2 points ! =6 match second derivative to first 3 points ! =7 match third derivative to first 4 points ! =else use knot-a-knot ! n-number of data points or knots-(n.ge.2) ! x(n)-abscissas of the knots in strictly increasing order ! f(1,i)-ordinates of the knots ! f(2,1)-input value of s'(x1) for k_bc1=1 ! f(2,n)-input value of s'(xn) for k_bcn=1 ! f(3,1)-input value of s''(x1) for k_bc1=2 ! f(3,n)-input value of s''(xn) for k_bcn=2 ! wk(n)-scratch work area for periodic BC !Output: ! f(2,i)=s'(x(i)) ! f(3,i)=s''(x(i)) ! f(4,i)=s'''(x(i))

Home | Top |

Hybrid

[Added to pspline -- DMC April 2007] The Pspline library now allows the definition of hybrid interpolating functions for data sets of dimensionality 2 or higher. In a hybrid interpolant, the interpolation method is controlled independently along each dimension of the data. In this context there are 4 interpolation methods that can be mixed: -1: zonal step function --no boundary condition --number of data points = dimension grid size - 1 0: piecewise linear --no boundary condition --number of data points = dimension grid size 1: C1 Akima Hermite --boundary condition: default, periodic, or control 1st derivative. --number of data points = dimension grid size 2: C2 Compact Spline --boundary condition: default, periodic, control 1st or 2nd derivative. --number of data points = dimension grid size These interpolation methods are specified separately for each dimension of the data to be interpolated. => For technical reasons in the implementation, there is one restriction: cubic interpolation methods (1) and (2) cannot be mixed. Note that a zonal step function can be specified along 1 or more of the data dimensions. But, if this is used, the size of the data along zonal step function dimensions must be ONE LESS than the size of the corresponding grid (as used e.g. in genxpkg), which specifies the zone boundaries. Although hybrid setup and evaluation routines could be provided for datasets of dimensionality greater than 3, this has not been done as of April 2007. As with all other pspline library routines, the floating point arguments of the routines named below are of precision REAL. To get routines with arguments of precision REAL*8, prepend the prefix "R8" to the subroutine name. Routines for generating Hybrid data sets: ========================================= for 2d data sets f(i,j) = f(x(i),y(j)): mkintrp2d (real*8 version: r8mkintrp2d) for 3d data sets f(i,j,k) = f(x(i),y(j),z(k)): mkintrp3d (real*8 version: r8mkintrp3d) Vectorized Evaluation Routines (recommended): ============================================= given a location vector, do lookup & evaluation: vecintrp2d vecintrp3d given a grid (different than the original grid), to map the function from the old grid to the new grid: gridintrp2d gridintrp3d if zone location data is already computed, to do just the vectorized spline evaluation: fvintrp2d fvintrp3d Scalar Evaluation Routines ========================== evintrp2d evintrp3d

Home | Top |

Note_on_derivatives

Derivatives can be requested along any dimension, but, the value returned will be ZERO if the derivative order exceeds the capability of the interpolation method along that dimension. Thus: any derivative along a zonal step function dimension results in ZERO. A 2nd or higher derivative along a piecewise linear dimension yields ZERO. A 2nd or higher derivative along a Hermite dimension yields ZERO.

Home | Top |

mkintrp2d

subroutine mkintrp2d(x,nx,y,ny,jspline, > f,icoeff,ixdim,iydim, > ibcxmin,bcxmin,ibcxmax,bcxmax, > ibcymin,bcymin,ibcymax,bcymax, > ier) C C setup bicubic spline, or bicubic hermite, or Hybrid linear/zonal with C 1d cubic Hermite or Spline interpolation in 1 dimension C implicit NONE 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 :: jspline(2) ! interpolation method control C (1) -- 1st dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline C (2) -- 2nd dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline C C Standard interpolation-- all jspline values match C e.g. jspline(1)=jspline(2)=2 for bicubic spline C Hybrid interpolation occurs when multiple jspline values are specified. C C RESTRICTION: jspline = (/ 1, 2 /) or (/ 2, 1 /) not allowed. I.e. C Hermite and Spline interpolation cannot be mixed in a hybrid interpolant. C This restriction exists because of technical issues in the C implementation (it could be removed in principle but the work to do C this has not been scheduled). C C coefficient buffer dimensioning: integer :: icoeff ! #coefficients per data point integer :: ixdim ! nx; nx-1 if jspline(1)==-1 integer :: iydim ! ny; ny-1 if jspline(2)==-1 C input/output: real f(icoeff,ixdim,iydim) ! data & spline coefficients C C input bdy condition data: C meaning of arguments as described in "mkbicubw.for" C *ignored* for dimensions with piecewise linear or zonal interpolation; C For Hermite interpolation only the values {-1,0,1} are accepted. C integer ibcxmin ! bc flag for x=xmin real bcxmin(iydim) ! bc data vs. y at x=xmin integer ibcxmax ! bc flag for x=xmax real bcxmax(iydim) ! bc data vs. y at x=xmax C integer ibcymin ! bc flag for y=ymin real bcymin(ixdim) ! bc data vs. x at y=ymin integer ibcymax ! bc flag for y=ymax real bcymax(ixdim) ! bc data vs. x at y=ymax C integer ier ! =0 on exit if there is no error. c ier -- completion code, 0 for normal

Home | Top |

mkintrp3d

subroutine mkintrp3d(x,nx,y,ny,z,nz,jspline, > f,icoeff,ixdim,iydim,izdim, > ibcxmin,bcxmin,ibcxmax,bcxmax, > ibcymin,bcymin,ibcymax,bcymax, > ibczmin,bczmin,ibczmax,bczmax, > ier) c c setup a tricubic spline, or tricubic Hermite, or hybrid linear/zonal 2d or c 1d with 1d or 2d cubic or Hermite spline interpolation c implicit NONE 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 :: jspline(3) ! interpolation method control C (1) -- 1st dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline C (2) -- 2nd dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline C (3) -- 3rd dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline C C Standard interpolation-- all jspline values match C e.g. jspline(1)=jspline(2)=jspline(3)=2 for tricubic spline C C Hybrid interpolation-- not all jspline values are the same. C C RESTRICTION: if any jspline(...) element has value 1, none can have C value 2. I.e. Spline and Hermite interpolation cannot currently be mixed. C This restriction exists because of technical issues in the C implementation (it could be removed in principle but the work to do C this has not been scheduled). c c coefficient buffer dimensions integer :: icoeff ! #coefficients per data point integer :: ixdim ! nx; nx-1 if jspline(1)==-1 integer :: iydim ! ny; ny-1 if jspline(2)==-1 integer :: izdim ! nz; nz-1 if jspline(3)==-1 c input/output: real f(icoeff,ixdim,iydim,izdim) ! data and spline coefficients c C boundary condition data C Contiguous storage is assumed-- 1st dimension size must match C actual use C 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 C real bcxmin(iydim,izdim),bcxmax(iydim,izdim) ! xmin & xmax BC data real bcymin(ixdim,izdim),bcymax(ixdim,izdim) ! ymin & ymax BC data real bczmin(ixdim,iydim),bczmax(ixdim,iydim) ! zmin & zmax BC data 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 not required for zonal or piecewise linear interpolation c for Hermite interpolation ibc* values from set {-1,0,1} are accepted. 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 -- df/dx BC from 1st divided difference c =6 -- d2f/dx2 BC from 2nd divided difference (parabolic fit) c =7 -- d3f/dx3 BC from 3rd divided difference (cubic fit) 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 c interpretation of ibcymin,bcymin,ibcymax,bcymax c is same as with ibcxmin,... c c interpretation of ibczmin,bczmin,ibczmax,bczmax c is same as with ibcxmin,... c c the explicit bdy condition arrays are referenced only if the c corresponding "ibc" flag values are set > 0. c integer ier ! exit code c ier -- completion code, 0 for normal

Home | Top |

vecintrp2d

subroutine vecintrp2d(ict,ivec,xvec,yvec,ivd,fval, > nx,xpkg,ny,ypkg,jspline,fspl,icoeff,ixdim,iydim, > iwarn,ier) c implicit NONE c c vectorized hybrid spline evaluation routine -- 2d c 1. call vectorized zone lookup routine c 2. call vectorized hybrid 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 c ict(1)=3 followed by ict(2:5) = 1 or 0 allow 3rd derivatives to c be selected: fxxx fxxy fxyy fyyy c c ict(1)=4 followed by ict(2:4) allow 4th derivatives to c be selected: fxxxy fxxyy fxyyy; fxxxx=fyyyy=0 c c ict(1)=5 followed by ict(2:4) allow 4th derivatives to c be selected: fxxxyy fxxyyy c c ict(1)=6 specifies 6th derivative: fxxxyyy (step fcn) c c in hybrid spline evaluation, any derivative along a dimension c with zonal interpolation (jspline(i)=-1) gives zero; c c piecewise linear and Hermite interpolation give zero if a 2nd or c higher derivative is sought along any dimension. 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) c integer :: jspline(2) ! interpolation method on each dim: c jspline(1) for x; jspline(2) for y c -1: zonal step fcn; 0: pcwise linear; 1: Hermite; 2: compact spline c integer icoeff ! #coeffs per data node integer ixdim ! x dim for fspl integer iydim ! y dim for fspl real fspl(icoeff,ixdim,iydim) ! hybrid spline coefficients c c ixdim=nx-1 for zonal step function along x dimension; o.w. expect ixdim=nx c iydim=ny-1 for zonal step function along y dimension; o.w. expect iydim=ny 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

Home | Top |

vecintrp3d

subroutine vecintrp3d(ict,ivec,xvec,yvec,zvec,ivd,fval, > nx,xpkg,ny,ypkg,nz,zpkg,jspline,fspl,icoeff,ixdim,iydim,izdim, > iwarn,ier) c implicit NONE c c vectorized hybrid spline evaluation routine -- 3d c 1. call vectorized zone lookup routine c 2. call vectorized hybrid 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 higher derivatives can be selected by setting ict(1)=3,-3,4,-4,... c see fvtricub.for comments. c c in hybrid spline evaluation, any derivative along a dimension c with zonal interpolation (jspline(i)=-1) gives zero; c c piecewise linear and Hermite interpolation give zero if a 2nd or c higher derivative is sought along any dimension. 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) c integer :: jspline(3) ! interpolation method on each dim: c jspline(1) for x; jspline(2) for y; jspline(3) for z c -1: zonal step fcn; 0: pcwise linear; 1: Hermite; 2: compact spline c integer icoeff ! #coeffs per data node integer ixdim ! x dim for fspl integer iydim ! y dim for fspl integer izdim ! z dim for fspl real fspl(icoeff,ixdim,iydim,izdim) ! hybrid spline coefficients c c ixdim=nx-1 for zonal step function along x dimension; o.w. expect ixdim=nx c iydim=ny-1 for zonal step function along y dimension; o.w. expect iydim=ny c izdim=nz-1 for zonal step function along z dimension; o.w. expect izdim=nz 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

Home | Top |

gridintrp2d

subroutine gridintrp2d( > x_newgrid,nx_new,y_newgrid,ny_new,f_new,if1, > nx,xpkg,ny,ypkg,jspline,fspl,icoeff,ixdim,iydim,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 jspline(2) ! interpolation type by dimension integer :: icoeff ! coefficients per data point integer :: ixdim,iydim ! fspl dimensions ! =nx,ny unless zonal step function interpolation is used real fspl(icoeff,ixdim,iydim) ! 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

Home | Top |

gridintrp3d

subroutine gridintrp3d( > x_newgrid,nx_new,y_newgrid,ny_new,z_newgrid,nz_new, > f_new,if1,if2, > nx,xpkg,ny,ypkg,nz,zpkg, > jspline,fspl,icoeff,ixdim,iydim,izdim, > 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 :: jspline(3) ! interpolation type, by dimension ! -1: zonal step fcn; 0: pclin; 1: Hermite; 2: Spline integer :: icoeff ! coefficients per data point integer :: ixdim,iydim,izdim ! fspl dimensions ! ixdim=nx unless zonal step fcn interpolation in x is used ! similar comment applies to iydim,izdim real fspl(icoeff,ixdim,iydim,izdim) ! 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

Home | Top |

fvintrp2d

C evaluate Hybrid function interpolation -- 2d fcn C NO ERROR CHECKING C subroutine fvintrp2d(ict,ivec,ivecd, > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi, > jspline,fin,icoeff,ixdim,iydim) C implicit NONE 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 integer :: jspline(2) ! interpolation type control integer :: icoeff ! coefficient dimension integer :: ixdim,iydim ! x & y dimensions real fin(icoeff,ixdim,iydim) ! data on which to interpolate... C real fval(ivecd,6) ! output returned C C for detailed description of fin, ict and fval see subroutine C evintrp2d 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 C

Home | Top |

fvintrp3d

C evaluate Hybrid function interpolation -- 3d fcn C NO ERROR CHECKING C subroutine fvintrp3d(ict,ivec,ivecd, > fval,ii,jj,kk,xparam,yparam,zparam, > hx,hxi,hy,hyi,hz,hzi, > jspline,fin,icoeff,ixdim,iydim,izdim) C implicit NONE C C use mkintrp3d 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 integer :: jspline(3) ! interpolation type control integer :: icoeff ! coefficient dimension integer :: ixdim,iydim,izdim ! x & y & z dimensions real fin(icoeff,ixdim,iydim,izdim) ! data on which to interpolate... C real fval(ivecd,10) ! output returned C C for detailed description of fin, ict and fval see subroutine evintrp3d C 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 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

Home | Top |

evintrp2d

subroutine evintrp2d(xget,yget,x,nx,y,ny,jspline, > f,icoeff,ixdim,iydim,ict,fval,ier) C implicit NONE C C evaluate a 2d Hybrid interpolant on a rectilinear grid C C this subroutine calls three subroutines: C vecin2d_argchk -- error check C herm2xy -- find cell containing (xget,yget) C fvintrp2d -- evaluate interpolant function and (optionally) derivatives C C input arguments: C ================ C real xget,yget ! target of this interpolation integer nx,ny ! grid sizes real x(nx) ! ordered x grid real y(ny) ! ordered y grid C integer :: jspline(2) ! interpolation method for each C dimension: jspline(1) for x; jspline(2) for y C =-1: zonal step function; =0: pc lin; =1: Hermite; =2: Spline C integer :: icoeff ! no. of coefficients per data point integer :: ixdim,iydim ! dimensioning: ! ixdim=nx-1 for zonal step in x; otherwise nx ! iydim=ny-1 for zonal step in y; otherwise ny real f(icoeff,ixdim,iydim) ! function data w/coefficients C integer ict(6) ! code specifying output desired C C Note on derivatives: for dimensions along which zonal step function C interpolation is done, ANY derivative request returns ZERO. C For dimensions along which piecewise linear or Hermite interpolation C are done, more than one differentiation returns ZERO! C C Derivative controls are the same as for the compact 2d spline evaluation C routine (evbicub): 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 fifth output (depends on ict(...) spec) C fval(6) receives the sixth 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

Home | Top |

evintrp3d

subroutine evintrp3d(xget,yget,zget,x,nx,y,ny,z,nz,jspline, > f,icoeff,ixdim,iydim,izdim,ict,fval,ier) C implicit NONE C C use mkintrp3d to set up spline coefficients... C C evaluate a 3d Hybrid interpolant on a rectilinear grid C C this subroutine calls two subroutines: C vecin3d_argchk -- error check C herm3xyz -- find cell containing (xget,yget,zget) C fvintrp3d -- evaluate the spline function (w/derivatives if req.) C C input arguments: C ================ C real xget,yget,zget ! target of this interpolation integer nx,ny,nz ! grid sizes real x(nx) ! ordered x grid real y(ny) ! ordered y grid real z(nz) ! ordered z grid C integer :: jspline(3) ! interpolation method for each C dimension: jspline(1) for x; jspline(2) for y; jspline(3) for z C =-1: zonal step function; =0: pc lin; =1: Hermite; =2: Spline C integer :: icoeff ! no. of coefficients per data point integer :: ixdim,iydim,izdim ! dimensioning: ! ixdim=nx-1 for zonal step in x; otherwise nx ! iydim=ny-1 for zonal step in y; otherwise ny ! izdim=nz-1 for zonal step in z; otherwise nz C real f(icoeff,ixdim,iydim,izdim) ! function data C integer ict(10) ! code specifying output desired C C Note on derivatives: for dimensions along which zonal step function C interpolation is done, ANY derivative request returns ZERO. C For dimensions along which piecewise linear or Hermite interpolation C are done, more than one differentiation returns ZERO! C C Derivative controls are the same as for the compact 3d spline evaluation C routine (evtricub): 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 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

Home | Top |

Test_programs

The PSPLINE software is distributed with a standalone program `pspltest' which demonstrates the spline and hermite software by using it to fit analytic expressions in 1d, 2d, and 3d, and comparing the fit to the original. The location of the test program will depend on site specific installation decisions. TRANSP sites can find a copy of `pspltest' in the TRANSP executables area. `pspltest' runs the single precision test. For the R8 routines there is a double precision test program, `pspltes8'. Another standalone program, `qk_pspline', demonstrates the evaluation of all available derivatives of splines, spot tests that expected continuity conditions apply across cell boundaries for all function values and derivatives, and that the R8 precision "compact" splines (ezspline) and non-compact splines match (they should to near machine precision). Another pair of test programs exist for exercising the "ezspline" interface. These are: ezspline_test_r4, ezspline_test_r8. All the programs come with reference outputs. Comparing run outputs with reference outputs takes some care. Some output quantities, e.g. cpu timings, are bound to differ from machine to machine. And, binary results will not be bit for bit identical, so, outputs of "residuals" will show small differences from machine to machine.

Home | Top |

About this document

This Document was created byhlptohtml