Pittsburgh Supercomputing Center (PSC)BINT4
SUBROUTINE BINT4(X,Y,NDATA,IBCL,IBCR,FBCL,FBCR,KNTOPT,T,BCOEF,N, Written by D. E. Amos, August, 1979. References SAND-78-1968 A Practical Guide to Splines by C. de Boor, Applied Mathematics Series 27, Springer, 1979. SIAM J. Numerical Analysis, 14, No. 3, June 1977, pp. 441-472. Abstract BINT4 computes the B representation (T,BCOEF,N,K) of a cubic spline (K=4) which interpolates data (X(I)),Y(I))), I=1,NDATA. Parameters IBCL, IBCR, FBCL, FBCR allow the specification of the spline first or second derivative at both X(1) and X(NDATA). When this data is not specified by the problem, it is common practice to use a natural spline by setting second derivatives at X(1) and X(NDATA) to zero (IBCL=IBCR=2,FBCL=FBCR=0.0). The spline is defined on T(4) .LE. X .LE. T(N+1) with (ordered) interior knots at X(I)) values where N=NDATA+2. The knots T(1), T(2), T(3) lie to the left of T(4)=X(1) and the knots T(N+2), T(N+3), T(N+4) lie to the right of T(N+1)=X(NDATA) in increasing order. If no extrapolation outside (X(1),X(NDATA)) is anticipated, the knots T(1)=T(2)=T(3)=T(4)=X(1) and T(N+2)=T(N+3)=T(N+4)= T(N+1)=X(NDATA) can be specified by KNTOPT=1. KNTOPT=2 selects a knot placement for T(1), T(2), T(3) to make the first 7 knots symmetric about T(4)=X(1) and similarly for T(N+2), T(N+3), T(N+4) about T(N+1)=X(NDATA). KNTOPT=3 allows the user to make his own selection, in increasing order, for T(1), T(2), T(3) to the left of X(1) and T(N+2), T(N+3), T(N+4) to the right of X(NDATA) in the work array W(1) through W(6). In any case, the interpolation on T(4) .LE. X .LE. T(N+1) by using function BVALU is unique for given boundary conditions. BINT4 calls BSPVD, BNFAC, BNSLV, R1MACH, XERROR Description of Arguments Input X - X vector of abscissae of length NDATA, distinct and in increasing order Y - Y vector of ordinates of length NDATA NDATA - number of data points, NDATA .GE. 2 IBCL - selection parameter for left boundary condition IBCL = 1 constrain the first derivative at X(1) to FBCL = 2 constrain the second derivative at X(1) to FBCL IBCR - selection parameter for right boundary condition IBCR = 1 constrain first derivative at X(NDATA) to FBCR IBCR = 2 constrain second derivative at X(NDATA) to FBCR FBCL - left boundary values governed by IBCL FBCR - right boundary values governed by IBCR KNTOPT - knot selection parameter KNTOPT = 1 sets knot multiplicity at T(4) and T(N+1) to 4 = 2 sets a symmetric placement of knots about T(4) and T(N+1) = 3 sets TNP)=WNP) and T(N+1+I)=w(3+I),I=1,3 where WNP),I=1,6 is supplied by the user W - work array of dimension at least 5*(NDATA+2) if KNTOPT=3, then W(1),W(2),W(3) are knot values to the left of X(1) and W(4),W(5),W(6) are knot values to the right of X(NDATA) in increasing order to be supplied by the user Output T - knot array of length N+4 BCOEF - B-spline coefficient array of length N N - number of coefficients, N=NDATA+2 K - order of spline, K=4 Error Conditions Improper input is a fatal error Singular system of equations is a fatal error
BINTK
SUBROUTINE BINTK(X,Y,T,N,K,BCOEF,Q,WORK) Written by Carl de Boor and modified by D. E. Amos References A Practical Guide to Splines by C. de Boor, Applied Mathematics Series 27, Springer, 1979. Abstract BINTK is the SPLINT routine of the reference. BINTK produces the B-spline coefficients, BCOEF, of the B-spline of order K with knots T(I), I=1,...,N+K, which takes on the value Y(I) at X(I), I=1,...,N. The spline or any of its derivatives can be evaluated by calls to BVALU. The I-th equation of the linear system A*BCOEF = B for the coefficients of the interpolant enforces interpolation at X(I)), I=1,...,N. Hence, B(I) = Y(I), all I, and A is a band matrix with 2K-1 bands if A is invertible. The matrix A is generated row by row and stored, diagonal by diagonal, in the rows of Q, with the main diagonal going into row K. The banded system is then solved by a call to BNFAC (which constructs the triangular factorization for A and stores it again in Q), followed by a call to BNSLV (which then obtains the solution BCOEF by substitution). BNFAC does no pivoting, since the total positivity of the matrix A makes this unnecessary. The linear system to be solved is (theoretically) invertible if and only if T(I) .LT. X(I)) .LT. T(I+K), all I. Equality is permitted on the left for I=1 and on the right for I=N when K knots are used at X(1) or X(N). Otherwise, violation of this condition is certain to lead to an error. BINTK calls BSPVN, BNFAC, BNSLV, XERROR Description of Arguments Input X - vector of length N containing data point abscissa in strictly increasing order. Y - corresponding vector of length N containing data point ordinates. T - knot vector of length N+K since T(1),..,T(K) .LE. X(1) and T(N+1),..,T(N+K) .GE. X(N), this leaves only N-K knots (not nec- essarily X(I)) values) interior to (X(1),X(N)) N - number of data points, N .GE. K K - order of the spline, K .GE. 1 Output BCOEF - a vector of length N containing the B-spline coefficients Q - a work vector of length (2*K-1)*N, containing the triangular factorization of the coefficient matrix of the linear system being solved. The coefficients for the interpolant of an additional data set (X(I)),YY(I)), I=1,...,N with the same abscissa can be obtained by loading YY into BCOEF and then executing call BNSLV(Q,2K-1,N,K-1,K-1,BCOEF) WORK - work vector of length 2*K Error Conditions Improper input is a fatal error Singular system of equations is a fatal error
BSPDOC
SUBROUTINE BSPDOC Written by D. E. Amos, May 1980 References 1. Computation with Splines and B-Splines by D. E. Amos, SAND78-1968, March, 1979. 2. Quadrature Subroutines for Splines and B-Splines by D. E. Amos, SAND79-1825, December, 1979. 3. A Practical Guide to Splines by C. de Boor, Applied Math Sci 27, Springer, N.Y., 1978. 4. On Calculating with B-Splines by C. de Boor, J. Approx. Theory,6,50-62(1972) 5. Package for Calculating with B-Splines by C. De Boor, SIAM J. Numer. anal.,14,441-472(1977). 6. Constrained Least Squares Curve Fitting to Discrete Data Using B-Splines - A User's Guide by R. J. Hanson, SAND78-1291, February,1979. 7. Monotone Piecewise Cubic Interpolation by F. N. Fritsch and R. E. Carlson, LLNL report UCRL-82453, January, 1979. Abstract BSPDOC is a non-executable, B-spline documentary routine. The narrative describes a B-spline and the routines necessary to manipulate B-splines at a fairly high level. The basic package described herein is that of reference 5 with names altered to prevent duplication and conflicts with routines from reference 3. The call lists used here are also different. Work vectors were added to ensure portability and proper execution in an overlay environ- ment. These work arrays can be used for other purposes except as noted in BSPVN. While most of the original routines in reference 5 were restricted to orders 20 or less, this restriction was removed from all routines except the quadrature routine BSQAD. (See the section below on differentiation and integration for details.) The subroutines referenced below are single precision routines. Corresponding double precision versions are also part of the package, and these are referenced by prefixing a D in front of the single precision name. For example, BVALU and DBVALU are the single and double precision versions for evaluating a B-spline or any of its deriva- tives in the B-representation. ****Description of B-Splines**** A collection of polynomials of fixed degree K-1 defined on a subdivision (X(I),X(I+1)), I=1,...,M-1 of (A,B) with X(1)=A, X(M)=B is called a B-spline of order K. If the spline has K-2 continuous derivatives on (A,B), then the B-spline is simply called a spline of order K. Each of the M-1 polynomial pieces has K coefficients, making a total of K(M-1) parameters. This B-spline and its derivatives have M-2 jumps at the subdivision points X(I), I=2,...,M-1. Continuity requirements at these subdivision points add constraints and reduce the number of free parameters. If a B-spline is continuous at each of the M-2 sub- division points, there are K(M-1)-(M-2) free parameters; if in addition the B-spline has continuous first derivatives, there are K(M-1)-2(M-2) free parameters, etc., until we get to a spline where we have K(M-1)-(K-1)(M-2) = M+K-2 free parameters. Thus, the principle is that increasing the continuity of derivatives decreases the number of free parameters and conversely. The points at which the polynomials are tied together by the continuity conditions are called knots. If two knots are allowed to come together at some X(I), then we say that we have a knot of multiplicity 2 there, and the knot values are the X(I) value. If we reverse the procedure of the first paragraph, we find that adding a knot to increase multiplicity increases the number of free parameters and, according to the principle above, we thereby introduce a discontinuity in what was the highest continuous derivative at that knot. Thus, the number of free parameters is N = NU+K-2 where NU is the sum of multipicities at the X(I) values with X(1) and X(M) of multiplicity 1 (NU = M if all knots are simple, i.e., for a spline, all knots have multiplicity 1.) Each knot can have a multiplicity of at most K. A B-spline is commonly written in the B-representation Y(X) = sum( A(I)*B(I,X), I=1 , N) to show the explicit dependence of the spline on the free parameters or coefficients A(I)=BCOEF(I) and basis functions B(I,X). These basis functions are themselves special B-splines which are zero except on (at most) K adjoining intervals where each B(I,X) is positive and, in most cases, hat or bell- shaped. In order for the nonzero part of B(1,X) to be a spline covering (X(1),X(2)), it is necessary to put K-1 knots to the left of A and similarly for B(N,X) to the right of B. Thus, the total number of knots for this representation is NU+2K-2 = N+K. These knots are carried in an array T(*) dimensioned by at least N+K. From the construction, A=T(K) and B=T(N+1) and the spline is defined on T(K).LE.X.LE.T(N+1). The nonzero part of each basis function lies in the Interval (T(I),T(I+K)). In many problems where extrapolation beyond A or B is not anticipated, it is common practice to set T(1)=T(2)=...=T(K)=A and T(N+1)=T(N+2)=...= T(N+K)=B. In summary, since T(K) and T(N+1) as well as interior knots can have multiplicity K, the number of free parameters N = sum of multiplicties - K. The fact that each B(I,X) function is nonzero over at most K intervals means that for a given X value, there are at most K nonzero terms of the sum. This leads to banded matrices in linear algebra problems, and references 3 and 6 take advantage of this in con- structing higher level routines to achieve speed and avoid ill-conditioning. ****Basic Routines**** The basic routines which most casual users will need are those concerned with direct evaluation of splines or B-splines. Since the B-representation, denoted by (T,BCOEF,N,K), is preferred because of numerical stability, the knots T(*), the B-spline coefficients BCOEF(*), the number of coefficients N, and the order K of the polynomial pieces (of degree K-1) are usually given. While the knot array runs from T(1) to T(N+K), the B-spline is normally defined on the interval T(K).LE.X.LE. T(N+1). To evaluate the B-spline or any of its derivatives on this interval, one can use Y = BVALU(T,BCOEF,N,K,ID,X,INBV,WORK) where ID is an integer for the ID-th derivative, 0.LE.ID.LE.K-1. ID=0 gives the zero-th derivative or B-spline value at X. If X.LT.T(K) or X.GT.T(N+1), whether by mistake or the result of round off accumulation in incrementing X, BVALU gives a diagnostic. INBV is an initialization parameter which is set to 1 on the first call. Distinct splines require distinct INBV parameters. WORK is a scratch vector of length at least 3*K. When more conventional communication is needed for publication, physical interpretation, etc., the B-spline coefficients can be converted to piecewise polynomial (PP) coefficients. Thus, the breakpoints (distinct knots) XI(*), the number of polynomial pieces LXI, and the (right) derivatives C(*,J) at each breakpoint XI(J) are needed to define the Taylor expansion to the right of XI(J) on each interval XI(J).LE. X.LT.XI(J+1), J=1,LXI where XI(1)=A and XI(LXI+1)=B. These are obtained from the (T,BCOEF,N,K) representation by CALL BSPPP(T,BCOEF,N,K,LDC,C,XI,LXI,WORK) where LDC.GE.K is the leading dimension of the matrix C and WORK is a scratch vector of length at least K*(N+3). Then the PP-representation (C,XI,LXI,K) of Y(X), denoted by Y(J,X) on each interval XI(J).LE.X.LT.XI(J+1), is Y(J,X) = sum( C(I,J)*((X-XI(J))**(I-1))/factorial(I-1), I=1,K) for J=1,...,LXI. One must view this conversion from the B- to the PP-representation with some skepticism because the conversion may lose significant digits when the B-spline varies in an almost discontinuous fashion. To evaluate the B-spline or any of its derivatives using the PP- representation, one uses Y = PPVAL(LDC,C,XI,LXI,K,ID,X,INPPV) where ID and INPPV have the same meaning and useage as ID and INBV in BVALU. To determine to what extent the conversion process loses digits, compute the relative error ABS((Y1-Y2)/Y2) over the X interval with Y1 from PPVAL and Y2 from BVALU. A major reason for considering PPVAL is that evaluation is much faster than that from BVALU. Recall that when multiple knots are encountered, jump type discontinuities in the B-spline or its derivatives occur at these knots, and we need to know that BVALU and PPVAL return right limiting values at these knots except at X=B where left limiting values are returned. These values are used for the Taylor expansions about left end points of breakpoint intervals. That is, the derivatives C(*,J) are right derivatives. Note also that a computed X value which, mathematically, would be a knot value may differ from the knot by a round off error. When this happens in evaluating a dis- continuous B-spline or some discontinuous derivative, the value at the knot and the value at X can be radically different. In this case, setting X to a T or XI value makes the computation precise. For left limiting values at knots other than X=B, see the prologues to BVALU and other routines. ****Interpolation**** BINTK is used to generate B-spline parameters (T,BCOEF,N,K) which will interpolate the data by calls to BVALU. A similar interpolation can also be done for cubic splines using BINT4 or the code in reference 7. If the PP-representation is given, one can evaluate this representation at an appropriate number of abscissas to create data then use BINTK or BINT4 to generate the B-representation. ****Differentiation and Integration**** Derivatives of B-splines are obtained from BVALU or PPVAL. Integrals are obtained from BSQAD using the B-representation (T,BCOEF,N,K) and PPQAD using the PP-representation (C,XI,LXI, K). More compicated integrals involving the product of a of a function F and some derivative of a B-spline can be evaluated with BFQAD or PFQAD using the B- or PP- represen- tations respectively. All quadrature routines, except for PPQAD, are limited in accuracy to 18 digits or working precision, whichever is smaller. PPQAD is limited to working precision only. In addition, the order K for BSQAD is limited to 20 or less. If orders greater than 20 are required, use BFQAD with F(X) = 1. ****Extrapolation**** Extrapolation outside the interval (A,B) can be accomplished easily by the PP-representation using PPVAL. However, caution should be exercised, especially when several knots are located at A or B or when the extrapolation is carried significantly beyond A or B. On the other hand, direct evaluation with BVALU outside A=T(K).LE.X.LE.T(N+1)=B produces an error message, and some manipulation of the knots and coefficients are needed to extrapolate with BVALU. This process is described in reference 6. ****Curve Fitting and Smoothing**** Unless one has many accurate data points, direct inter- polation is not recommended for summarizing data. The results are often not in accordance with intuition since the fitted curve tends to oscillate through the set of points. Monotone splines (reference 7) can help curb this undulating tendency but constrained least squares is more likely to give an acceptable fit with fewer parameters. Subroutine FC, des- cribed in reference 6, is recommended for this purpose. The output from this fitting process is the B-representation. **** Routines in the B-Spline Package **** Single Precision Routines The subroutines referenced below are SINGLE PRECISION routines. Corresponding DOUBLE PRECISION versions are also part of the package and these are referenced by prefixing a D in front of the single precision name. For example, BVALU and DBVALU are the SINGLE and DOUBLE PRECISION versions for evaluating a B-spline or any of its deriva- tives in the B-representation. BINT4 - interpolates with splines of order 4 BINTK - interpolates with splines of order k BSQAD - integrates the B-representation on subintervals PPQAD - integrates the PP-representation BFQAD - integrates the product of a function F and any spline derivative in the B-representation PFQAD - integrates the product of a function F and any spline derivative in the PP-representation BVALU - evaluates the B-representation or a derivative PPVAL - evaluates the PP-representation or a derivative INTRV - gets the largest index of the knot to the left of x BSPPP - converts from B- to PP-representation BSPVD - computes nonzero basis functions and derivatives at x BSPDR - sets up difference array for BSPEV BSPEV - evaluates the B-representation and derivatives BSPVN - called by BSPEV, BSPVD, BSPPP and BINTK for function and derivative evaluations Auxiliary Routines BSGQ8,PPGQ8,BNSLV,BNFAC,XERROR,DBSGQ8,DPPGQ8,DBNSLV,DBNFAC Machine Dependent Routines I1MACH, R1MACH, D1MACH
BSPDR
SUBROUTINE BSPDR(T,A,N,K,NDERIV,AD) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract BSPDR is the BSPLDR routine of the reference. BSPDR uses the B-representation (T,A,N,K) to construct a divided difference table ADIF preparatory to a (right) derivative calculation in BSPEV. The lower triangular matrix ADIF is stored in vector AD by columns. The arrays are related by ADIF(I,J) = AD(I-J+1 + (2*N-J+2)*(J-1)/2) I = J,N , J = 1,NDERIV . Description of Arguments Input T - knot vector of length N+K A - B-spline coefficient vector of length N N - number of B-spline coefficients N = sum of knot multiplicities-K K - order of the spline, K .GE. 1 NDERIV - number of derivatives, 1 .LE. NDERIV .LE. K. NDERIV=1 gives the zero-th derivative = function value Output AD - table of differences in a vector of length (2*N-NDERIV+1)*NDERIV/2 for input to BSPEV Error Conditions Improper input is a fatal error
BSPEV
SUBROUTINE BSPEV(T,AD,N,K,NDERIV,X,INEV,SVALUE,WORK) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract BSPEV is the BSPLEV routine of the reference. BSPEV calculates the value of the spline and its derivatives at X from the B-representation (T,A,N,K) and returns them in SVALUE(I),I=1,NDERIV, T(K) .LE. X .LE. T(N+1). AD(I) can be the B-spline coefficients A(I), I=1,N if NDERIV=1. Other- wise AD must be computed before hand by a call to BSPDR (T,A, N,K,NDERIV,AD). If X=T(I),I=K,N, right limiting values are obtained. To compute left derivatives or left limiting values at a knot T(I), replace N by I-1 and set X=T(I), I=K+1,N+1. BSPEV calls INTRV, BSPVN Description of Arguments Input T - knot vector of length N+K AD - vector of length (2*N-NDERIV+1)*NDERIV/2 containing the difference table from BSPDR. N - number of B-spline coefficients N = sum of knot multiplicities-K K - order of the B-spline, K .GE. 1 NDERIV - number of derivatives, 1 .LE. NDERIV .LE. K. NDERIV=1 gives the zero-th derivative = function value X - argument, T(K) .LE. X .LE. T(N+1) INEV - an initialization parameter which must be set to 1 the first time BSPEV is called. Output INEV - INEV contains information for efficient process- ing after the initial call and INEV must not be changed by the user. Distinct splines require distinct INEV parameters. SVALUE - vector of length NDERIV containing the spline value in SVALUE(1) and the NDERIV-1 derivatives in the remaining components. WORK - work vector of length 3*K Error Conditions Improper input is a fatal error.
BSPLVD
SUBROUTINE BSPLVD(T,K,X,ILEFT,VNIKX,NDERIV)
BSPLVN
SUBROUTINE BSPLVN(T,JHIGH,INDEX,X,ILEFT,VNIKX)
BSPPP
SUBROUTINE BSPPP(T,A,N,K,LDC,C,XI,LXI,WORK) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract BSPPP is the BSPLPP routine of the reference. BSPPP converts the B-representation (T,A,N,K) to the piecewise polynomial (PP) form (C,XI,LXI,K) for use with PPVAL. Here XI(*), the break point array of length LXI, is the knot array T(*) with multiplicities removed. The columns of the matrix C(I,J) contain the right Taylor derivatives for the polynomial expansion about XI(J) for the intervals XI(J) .LE. X .LE. XI(J+1), I=1,K, J=1,LXI. Function PPVAL makes this evaluation at a specified point X in XI(1) .LE. X .LE. XI(LXI(1) .LE. X .LE. XI+1) BSPPP calls BSPDR, BSPEV, INTRV, BSPVN Description of Arguments Input T - knot vector of length N+K A - B-spline coefficient vector of length N N - number of B-spline coefficients N = sum of knot multiplicities-K K - order of the B-spline, K .GE. 1 LDC - leading dimension of C, LDC .GE. K Output C - matrix of dimension at least (K,LXI) containing right derivatives at break points XI - XI break point vector of length LXI+1 LXI - number of break points, LXI .LE. N-K+1 WORK - work vector of length K*(N+3) Error Conditions Improper input is a fatal error
BSPVD
SUBROUTINE BSPVD(T,K,NDERIV,X,ILEFT,LDVNIK,VNIKX,WORK) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract BSPVD is the BSPLVD routine of the reference. BSPVD calculates the value and all derivatives of order less than NDERIV of all basis functions which do not (possibly) vanish at X. ILEFT is input such that T(ILEFT) .LE. X .LT. T(ILEFT+1). A call to INTRV(T,N+1,X, ILO,ILEFT,MFLAG) will produce the proper ILEFT. The output of BSPVD is a matrix VNIKX(I,J) of dimension at least (K,NDERIV) whose columns contain the K nonzero basis functions and their NDERIV-1 right derivatives at X, I=1,K, J=1,NDERIV. These basis functions have indices ILEFT-K+I, I=1,K, K .LE. ILEFT .LE. N. The nonzero part of the I-th basis function lies in (T(I),T(I+K)), I=1,N. If X=T(ILEFT+1) then VNIKX contains left limiting values (left derivatives) at T(ILEFT+1). In particular, ILEFT = N produces left limiting values at the right end point X=T(N+1). To obtain left limiting values at T(I), I=K+1,N+1, set X= next lower distinct knot, call INTRV to get ILEFT, set X=T(I), and then call BSPVD. BSPVD calls BSPVN Description of Arguments Input T - knot vector of length N+K, where N = number of B-spline basis functions N = sum of knot multiplicities-K K - order of the B-spline, K .GE. 1 NDERIV - number of derivatives = NDERIV-1, 1 .LE. NDERIV .LE. K X - argument of basis functions, T(K) .LE. X .LE. T(N+1) ILEFT - largest integer such that T(ILEFT) .LE. X .LT. T(ILEFT+1) LDVNIK - leading dimension of matrix VNIKX Output VNIKX - matrix of dimension at least (K,NDERIV) contain- ing the nonzero basis functions at X and their derivatives columnwise. WORK - a work vector of length (K+1)*(K+2)/2 Error Conditions Improper input is a fatal error
BSPVN
SUBROUTINE BSPVN(T,JHIGH,K,INDEX,X,ILEFT,VNIKX,WORK,IWORK) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract BSPVN is the BSPLVN routine of the reference. BSPVN calculates the value of all (possibly) nonzero basis functions at X of order MAX(JHIGH,(J+1)*(INDEX-1)), where T(K) .LE. X .LE. T(N+1) and J=IWORK is set inside the routine on the first call when INDEX=1. ILEFT is such that T(ILEFT) .LE. X .LT. T(ILEFT+1). A call to INTRV(T,N+1,X,ILO,ILEFT, MFLAG) produces the proper ILEFT. BSPVN calculates using the basic algorithm needed in BSPVD. If only basis functions are desired, setting JHIGH=K and INDEX=1 can be faster than calling BSPVD, but extra coding is required for derivatives (INDEX=2) and BSPVD is set up for this purpose. Left limiting values are set up as described in BSPVD. Description of Arguments Input T - knot vector of length N+K, where N = number of B-spline basis functions N = sum of knot multiplicities-K JHIGH - order of B-spline, 1 .LE. JHIGH .LE. K K - highest possible order INDEX - INDEX = 1 gives basis functions of order JHIGH = 2 denotes previous entry with WORK, IWORK values saved for subsequent calls to BSPVN. X - argument of basis functions, T(K) .LE. X .LE. T(N+1) ILEFT - largest integer such that T(ILEFT) .LE. X .LT. T(ILEFT+1) Output VNIKX - vector of length K for spline values. WORK - a work vector of length 2*K IWORK - a work parameter. Both WORK and IWORK contain information necessary to continue for INDEX = 2. When INDEX = 1 exclusively, these are scratch variables and can be used for other purposes. Error Conditions Improper input is a fatal error.
BVALU
FUNCTION BVALU(T,A,N,K,IDERIV,X,INBV,WORK) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract BVALU is the BVALUE function of the reference. BVALU evaluates the B-representation (T,A,N,K) of a B-spline at X for the function value on IDERIV = 0 or any of its derivatives on IDERIV = 1,2,...,K-1. Right limiting values (right derivatives) are returned except at the right end point X=T(N+1) where left limiting values are computed. The spline is defined on T(K) .LE. X .LE. T(N+1). BVALU returns a fatal error message when X is outside of this interval. To compute left derivatives or left limiting values at a knot T(I), replace N by I-1 and set X=T(I), I=K+1,N+1. BVALU calls INTRV Description of Arguments Input T - knot vector of length N+K A - B-spline coefficient vector of length N N - number of B-spline coefficients N = sum of knot multiplicities-K K - order of the B-spline, K .GE. 1 IDERIV - order of the derivative, 0 .LE. IDERIV .LE. K-1 IDERIV=0 returns the B-spline value X - argument, T(K) .LE. X .LE. T(N+1) INBV - an initialization parameter which must be set to 1 the first time BVALU is called. Output INBV - INBV contains information for efficient process- ing after the initial call and INBV must not be changed by the user. Distinct splines require distinct INBV parameters. WORK - work vector of length 3*K. BVALU - value of the IDERIV-th derivative at X Error Conditions An improper input is a fatal error
BVALUE
FUNCTION BVALUE(T,A,N,K,X,IDERIV)
BVDER
SUBROUTINE BVDER(X,Y,YP,G,IPAR)
BVPOR
SUBROUTINE BVPOR(Y,NROWY,NCOMP,XPTS,NXPTS,A,NROWA,ALPHA,NIC,B,
CHFDV
SUBROUTINE CHFDV(X1,X2,F1,F2,D1,D2,NE,XE,FE,DE,NEXT,IERR) CHFDV: Cubic Hermite Function and Derivative Evaluator Evaluates the cubic polynomial determined by function values F1,F2 and derivatives D1,D2 on interval (X1,X2), together with its first derivative, at the points XE(J), J=1(1)NE. If only function values are required, use CHFEV, instead.
Calling sequence: INTEGER NE, NEXT(2), IERR REAL X1, X2, F1, F2, D1, D2, XE(NE), FE(NE), DE(NE) CALL CHFDV (X1,X2, F1,F2, D1,D2, NE, XE, FE, DE, NEXT, IERR) Parameters: X1,X2 -- (input) endpoints of interval of definition of cubic. (Error return if X1.EQ.X2 .) F1,F2 -- (input) values of function at X1 and X2, respectively. D1,D2 -- (input) values of derivative at X1 and X2, respectively. NE -- (input) number of evaluation points. (Error return if NE.LT.1 .) XE -- (input) real array of points at which the functions are to be evaluated. If any of the XE are outside the interval [X1,X2], a warning error is returned in NEXT. FE -- (output) real array of values of the cubic function defined by X1,X2, F1,F2, D1,D2 at the points XE. DE -- (output) real array of values of the first derivative of the same function at the points XE. NEXT -- (output) integer array indicating number of extrapolation points: NEXT(1) = number of evaluation points to left of interval. NEXT(2) = number of evaluation points to right of interval. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable" errors: IERR = -1 if NE.LT.1 . IERR = -2 if X1.EQ.X2 . (Output arrays have not been changed in either case.)CHFEV
SUBROUTINE CHFEV(X1,X2,F1,F2,D1,D2,NE,XE,FE,NEXT,IERR) CHFEV: Cubic Hermite Function EValuator Evaluates the cubic polynomial determined by function values F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points XE(J), J=1(1)NE.
Calling sequence: INTEGER NE, NEXT(2), IERR REAL X1, X2, F1, F2, D1, D2, XE(NE), FE(NE) CALL CHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR) Parameters: X1,X2 -- (input) endpoints of interval of definition of cubic. (Error return if X1.EQ.X2 .) F1,F2 -- (input) values of function at X1 and X2, respectively. D1,D2 -- (input) values of derivative at X1 and X2, respectively. NE -- (input) number of evaluation points. (Error return if NE.LT.1 .) XE -- (input) real array of points at which the function is to be evaluated. If any of the XE are outside the interval [X1,X2], a warning error is returned in NEXT. FE -- (output) real array of values of the cubic function defined by X1,X2, F1,F2, D1,D2 at the points XE. NEXT -- (output) integer array indicating number of extrapolation points: NEXT(1) = number of evaluation points to left of interval. NEXT(2) = number of evaluation points to right of interval. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable" errors: IERR = -1 if NE.LT.1 . IERR = -2 if X1.EQ.X2 . (The FE-array has not been changed in either case.)CHFIV
REAL FUNCTION CHFIV(X1,X2,F1,F2,D1,D2,A,B,IERR) CHFIV: Cubic Hermite Function Integral Evaluator. Called by PCHIA to evaluate the integral of a single cubic (in Hermite form) over an arbitrary interval (A,B).
Calling sequence: INTEGER IERR REAL X1, X2, F1, F2, D1, D2, A, B REAL VALUE, CHFIV VALUE = CHFIV (X1, X2, F1, F2, D1, D2, A, B, IERR) Parameters: VALUE -- (output) VALUE of the requested integral. X1,X2 -- (input) endpoints if interval of definition of cubic. (Must be distinct. Error return if not.) F1,F2 -- (input) function values at the ends of the interval. D1,D2 -- (input) derivative values at the ends of the interval. A,B -- (input) endpoints of interval of integration. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable errors": IERR = -1 if X1.EQ.X2 . (VALUE has not been set in this case.)CHFMC
INTEGER FUNCTION CHFMC(D1,D2,DELTA) CHFMC: Cubic Hermite Function Monotonicity Checker. Called by PCHMC to determine the monotonicity properties of the cubic with boundary derivative values D1,D2 and chord slope DELTA.
Calling sequence: DOUBLE PRECISION D1, D2, DELTA INTEGER ISMON ISMON = CHFMC (D1, D2, DELTA) Parameters: D1,D2 -- (input) derivative values at the ends of an interval. DELTA -- (input) data slope over that interval. ISMON -- (output) integer function value, indicating the monoto- nicity of the cubic segment: ISMON = -1 if function is strictly decreasing; ISMON = 0 if function is constant; ISMON = 1 if function is strictly increasing; ISMON = 2 if function is non-monotonic; ISMON = 3 if unable to determine. Fortran intrinsics used: SIGN. Other routines used: R1MACH.DBINT4
SUBROUTINE DBINT4(X,Y,NDATA,IBCL,IBCR,FBCL,FBCR,KNTOPT,T,BCOEF,N, Written by D. E. Amos, August, 1979. References SAND-78-1968 A Practical Guide to Splines by C. de Boor, Applied Mathematics Series 27, Springer, 1979. SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract **** a double precision routine **** DBINT4 computes the B representation (T,BCOEF,N,K) of a cubic spline (K=4) which interpolates data (X(I),Y(I)), I=1,NDATA. Parameters IBCL, IBCR, FBCL, FBCR allow the specification of the spline first or second derivative at both X(1) and X(NDATA). When this data is not specified by the problem, it is common practice to use a natural spline by setting second derivatives at X(1) and X(NDATA) to zero (IBCL=IBCR=2,FBCL=FBCR=0.0). The spline is defined on T(4) .LE. X .LE. T(N+1) with (ordered) interior knots at X(I) values where N=NDATA+2. The knots T(1),T(2),T(3) lie to the left of T(4)=X(1) and the knots T(N+2), T(N+3), T(N+4) lie to the right of T(N+1)=X(NDATA) in increasing order. If no extrapolation outside (X(1),X(NDATA)) is anticipated, the knots T(1)=T(2)=T(3)=T(4)=X(1) and T(N+2)=T(N+3)=T(N+4)= T(N+1)=X(NDATA) can be specified by KNTOPT=1. KNTOPT=2 selects a knot placement for T(1), T(2), T(3) to make the first 7 knots symmetric about T(4)=X(1) and similarly for T(N+2), T(N+3), T(N+4) about T(N+1)=X(NDATA). KNTOPT=3 allows the user to make his own selection, in increasing order, for T(1), T(2), T(3) to the left of X(1) and T(N+2), T(N+3), T(N+4) to the right of X(NDATA) in the work array W(1) through W(6). In any case, the interpolation on T(4) .LE. X .LE. T(N+1) by using function DBVALU is unique for given boundary conditions. DBINT4 calls DBSPVD, DBNFAC, DBNSLV, D1MACH, XERROR Description of Arguments Input X,Y,FBCL,FBCR,W are double precision X - X vector of abscissae of length NDATA, distinct and in increasing order Y - Y vector of ordinates of length NDATA NDATA - number of data points, NDATA .GE. 2 IBCL - selection parameter for left boundary condition IBCL = 1 constrain the first derivative at X(1) to FBCL = 2 constrain the second derivative at X(1) to FBCL IBCR - selection parameter for right boundary condition IBCR = 1 constrain first derivative at X(NDATA) to FBCR IBCR = 2 constrain second derivative at X(NDATA) to FBCR FBCL - left boundary values governed by IBCL FBCR - right boundary values governed by IBCR KNTOPT - knot selection parameter KNTOPT = 1 sets knot multiplicity at T(4) and T(N+1) to 4 = 2 sets a symmetric placement of knots about T(4) and T(N+1) = 3 sets T(I)=W(I) and T(N+1+I)=W(3+I),I=1,3 where W(I),I=1,6 is supplied by the user W - work array of dimension at least 5*(NDATA+2) If KNTOPT=3, then W(1),W(2),W(3) are knot values to the left of X(1) and W(4),W(5),W(6) are knot values to the right of X(NDATA) in increasing order to be supplied by the user Output T,BCOEF are double precision T - knot array of length N+4 BCOEF - B spline coefficient array of length N N - number of coefficients, N=NDATA+2 K - order of spline, K=4 Error Conditions Improper input is a fatal error Singular system of equations is a fatal error
DBINTK
SUBROUTINE DBINTK(X,Y,T,N,K,BCOEF,Q,WORK) Written by Carl de Boor and modified by D. E. Amos References A Practical Guide to Splines by C. de Boor, Applied Mathematics Series 27, Springer, 1979. Abstract **** a double precision routine **** DBINTK is the SPLINT routine of the reference. DBINTK produces the B-spline coefficients, BCOEF, of the B-spline of order K with knots T(I), I=1,...,N+K, which takes on the value Y(I) at X(I), I=1,...,N. The spline or any of its derivatives can be evaluated by calls to DBVALU. The I-th equation of the linear system A*BCOEF = B for the coefficients of the interpolant enforces interpolation at X(I), I=1,...,N. Hence, B(I) = Y(I), for all I, and A is a band matrix with 2K-1 bands if A is invertible. The matrix A is generated row by row and stored, diagonal by diagonal, in the rows of Q, with the main diagonal going into row K. The banded system is then solved by a call to DBNFAC (which constructs the triangular factorization for A and stores it again in Q), followed by a call to DBNSLV (which then obtains the solution BCOEF by substitution). DBNFAC does no pivoting, since the total positivity of the matrix A makes this unnecessary. The linear system to be solved is (theoretically) invertible if and only if T(I) .LT. X(I) .LT. T(I+K), for all I. Equality is permitted on the left for I=1 and on the right for I=N when K knots are used at X(1) or X(N). Otherwise, violation of this condition is certain to lead to an error. DBINTK calls DBSPVN, DBNFAC, DBNSLV, XERROR Description of Arguments Input X,Y,T are double precision X - vector of length N containing data point abscissa in strictly increasing order. Y - corresponding vector of length N containing data point ordinates. T - knot vector of length N+K Since T(1),..,T(K) .LE. X(1) and T(N+1),..,T(N+K) .GE. X(N), this leaves only N-K knots (not nec- essarily X(I) values) interior to (X(1),X(N)) N - number of data points, N .GE. K K - order of the spline, K .GE. 1 Output BCOEF,Q,WORK are double precision BCOEF - a vector of length N containing the B-spline coefficients Q - a work vector of length (2*K-1)*N, containing the triangular factorization of the coefficient matrix of the linear system being solved. The coefficients for the interpolant of an additional data set (X(I),yY(I)), I=1,...,N with the same abscissa can be obtained by loading YY into BCOEF and then executing CALL DBNSLV(Q,2K-1,N,K-1,K-1,BCOEF) WORK - work vector of length 2*K Error Conditions Improper input is a fatal error Singular system of equations is a fatal error
DBKIAS
SUBROUTINE DBKIAS(X,N,KTRMS,T,ANS,IND,MS,GMRN,H,IERR) DBKIAS computes repeated integrals of the K0 Bessel function by the asymptotic expansion
DBKISR
SUBROUTINE DBKISR(X,N,SUM,IERR) DBKISR computes repeated integrals of the K0 Bessel function by the series for N=0,1, and 2.
DBKSOL
SUBROUTINE DBKSOL(N,A,X)
DBSPDR
SUBROUTINE DBSPDR(T,A,N,K,NDERIV,AD) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract **** a double precision routine **** DBSPDR is the BSPLDR routine of the reference. DBSPDR uses the B-representation (T,A,N,K) to construct a divided difference table ADIF preparatory to a (right) derivative calculation in DBSPEV. The lower triangular matrix ADIF is stored in vector AD by columns. The arrays are related by ADIF(I,J) = AD(I-J+1 + (2*N-J+2)*(J-1)/2) I = J,N , J=1,NDERIV. Description of Arguments Input T,A are double precision T - knot vector of length N+K A - B-spline coefficient vector of length N N - number of B-spline coefficients N = sum of knot multiplicities-K K - order of the spline, K .GE. 1 NDERIV - number of derivatives, 1 .LE. NDERIV .LE. K. NDERIV=1 gives the zero-th derivative = function value Output AD is double precision AD - table of differences in a vector of length (2*N-NDERIV+1)*NDERIV/2 for input to DBSPEV Error Conditions Improper input is a fatal error
DBSPEV
SUBROUTINE DBSPEV(T,AD,N,K,NDERIV,X,INEV,SVALUE,WORK) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract **** a double precision routine **** DBSPEV is the BSPLEV routine of the reference. DBSPEV calculates the value of the spline and its derivatives at X from the B-representation (T,A,N,K) and returns them in SVALUE(I),I=1,NDERIV, T(K) .LE. X .LE. T(N+1). AD(I) can be the B-spline coefficients A(I), I=1,N) if NDERIV=1. Otherwise AD must be computed before hand by a call to DBSPDR (T,A,N,K, NDERIV,AD). If X=T(I),I=K,N), right limiting values are obtained. To compute left derivatives or left limiting values at a knot T(I), replace N by I-1 and set X=T(I), I=K+1,N+1. DBSPEV calls DINTRV, DBSPVN Description of Arguments Input T,AD,X, are double precision T - knot vector of length N+K AD - vector of length (2*N-NDERIV+1)*NDERIV/2 containing the difference table from DBSPDR. N - number of B-spline coefficients N = sum of knot multiplicities-K K - order of the B-spline, K .GE. 1 NDERIV - number of derivatives, 1 .LE. NDERIV .LE. K. NDERIV=1 gives the zero-th derivative = function value X - argument, T(K) .LE. X .LE. T(N+1) INEV - an initialization parameter which must be set to 1 the first time DBSPEV is called. Output SVALUE,WORK are double precision INEV - INEV contains information for efficient process- ing after the initial call and INEV must not be changed by the user. Distinct splines require distinct INEV parameters. SVALUE - vector of length NDERIV containing the spline value in SVALUE(1) and the NDERIV-1 derivatives in the remaining components. WORK - work vector of length 3*K Error Conditions Improper input is a fatal error.
DBSPPP
SUBROUTINE DBSPPP(T,A,N,K,LDC,C,XI,LXI,WORK) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract **** a double precision routine **** DBSPPP is the BSPLPP routine of the reference. DBSPPP converts the B-representation (T,A,N,K) to the piecewise polynomial (PP) form (C,XI,LXI,K) for use with DPPVAL. Here XI(*), the break point array of length LXI, is the knot array T(*) with multiplicities removed. The columns of the matrix C(I,J) contain the right Taylor derivatives for the polynomial expansion about XI(J) for the intervals XI(J) .LE. X .LE. XI(J+1), I=1,K, J=1,LXI. Function DPPVAL makes this evaluation at a specified point X in XI(1) .LE. X .LE. XI(LXI+1) DBSPPP calls DBSPDR, DBSPEV, DINTRV, DBSPVN Description of Arguments Input T,A are double precision T - knot vector of length N+K A - B-spline coefficient vector of length N N - number of B-spline coefficients N = sum of knot multiplicities-K K - order of the B-spline, K .GE. 1 LDC - leading dimension of C, LDC .GE. K Output C,XI,WORK are double precision C - matrix of dimension at least (K,LXI) containing right derivatives at break points XI - XI break point vector of length LXI+1 LXI - number of break points, LXI .LE. N-K+1 WORK - work vector of length K*(N+3) Error Conditions Improper input is a fatal error
DBSPVD
SUBROUTINE DBSPVD(T,K,NDERIV,X,ILEFT,LDVNIK,VNIKX,WORK) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract **** a double precision routine **** DBSPVD is the BSPLVD routine of the reference. DBSPVD calculates the value and all derivatives of order less than NDERIV of all basis functions which do not (possibly) vanish at X. ILEFT is input such that T(ILEFT) .LE. X .LT. T(ILEFT+1). A call to INTRV(T,N+1,X, ILO,ILEFT,MFLAG) will produce the proper ILEFT. The output of DBSPVD is a matrix VNIKX(I,J) of dimension at least (K,NDERIV) whose columns contain the K nonzero basis functions and their NDERIV-1 right derivatives at X, I=1,K, J=1,NDERIV. These basis functions have indices ILEFT-K+I, I=1,K, K .LE. ILEFT .LE. N. The nonzero part of the I-th basis function lies in (T(I),T(I+K)), I=1,N). If X=T(ILEFT+1) then VNIKX contains left limiting values (left derivatives) at T(ILEFT+1). In particular, ILEFT = N produces left limiting values at the right end point X=T(N+1). To obtain left limiting values at T(I), I=K+1,N+1, set X= next lower distinct knot, call INTRV to get ILEFT, set X=T(I), and then call DBSPVD. DBSPVD calls DBSPVN Description of Arguments Input T,X are double precision T - knot vector of length N+K, where N = number of B-spline basis functions N = sum of knot multiplicities-K K - order of the B-spline, K .GE. 1 NDERIV - number of derivatives = NDERIV-1, 1 .LE. NDERIV .LE. K X - argument of basis functions, T(K) .LE. X .LE. T(N+1) ILEFT - largest integer such that T(ILEFT) .LE. X .LT. T(ILEFT+1) LDVNIK - leading dimension of matrix VNIKX Output VNIKX,WORK are double precision VNIKX - matrix of dimension at least (K,NDERIV) contain- ing the nonzero basis functions at X and their derivatives columnwise. WORK - a work vector of length (K+1)*(K+2)/2 Error Conditions Improper input is a fatal error
DBSPVN
SUBROUTINE DBSPVN(T,JHIGH,K,INDEX,X,ILEFT,VNIKX,WORK,IWORK) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract **** a double precision routine **** DBSPVN is the BSPLVN routine of the reference. DBSPVN calculates the value of all (possibly) nonzero basis functions at X of order MAX(JHIGH,(J+1)*(INDEX-1)), where T(K) .LE. X .LE. T(N+1) and J=IWORK is set inside the routine on the first call when INDEX=1. ILEFT is such that T(ILEFT) .LE. X .LT. T(ILEFT+1). A call to DINTRV(T,N+1,X,ILO,ILEFT,MFLAG) produces the proper ILEFT. DBSPVN calculates using the basic algorithm needed in DBSPVD. If only basis functions are desired, setting JHIGH=K and INDEX=1 can be faster than calling DBSPVD, but extra coding is required for derivatives (INDEX=2) and DBSPVD is set up for this purpose. Left limiting values are set up as described in DBSPVD. Description of Arguments Input T,X are double precision T - knot vector of length N+K, where N = number of B-spline basis functions N = sum of knot multiplicities-K JHIGH - order of B-spline, 1 .LE. JHIGH .LE. K K - highest possible order INDEX - INDEX = 1 gives basis functions of order JHIGH = 2 denotes previous entry with work, IWORK values saved for subsequent calls to DBSPVN. X - argument of basis functions, T(K) .LE. X .LE. T(N+1) ILEFT - largest integer such that T(ILEFT) .LE. X .LT. T(ILEFT+1) Output VNIKX, WORK are double precision VNIKX - vector of length K for spline values. WORK - a work vector of length 2*K IWORK - a work parameter. Both WORK and IWORK contain information necessary to continue for INDEX = 2. When INDEX = 1 exclusively, these are scratch variables and can be used for other purposes. Error Conditions Improper input is a fatal error.
DBVALU
DOUBLE PRECISION FUNCTION DBVALU(T,A,N,K,IDERIV,X,INBV,WORK) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract **** a double precision routine **** DBVALU is the BVALUE function of the reference. DBVALU evaluates the B-representation (T,A,N,K) of a B-spline at X for the function value on IDERIV=0 or any of its derivatives on IDERIV=1,2,...,K-1. Right limiting values (right derivatives) are returned except at the right end point X=T(N+1) where left limiting values are computed. The spline is defined on T(K) .LE. X .LE. T(N+1). DBVALU returns a fatal error message when X is outside of this interval. To compute left derivatives or left limiting values at a knot T(I), replace N by I-1 and set X=T(I), I=K+1,N+1. DBVALU calls DINTRV Description of Arguments Input T,A,X are double precision T - knot vector of length N+K A - B-spline coefficient vector of length N N - number of B-spline coefficients N = sum of knot multiplicities-K K - order of the B-spline, K .GE. 1 IDERIV - order of the derivative, 0 .LE. IDERIV .LE. K-1 IDERIV = 0 returns the B-spline value X - argument, T(K) .LE. X .LE. T(N+1) INBV - an initialization parameter which must be set to 1 the first time DBVALU is called. Output WORK,DBVALU are double precision INBV - INBV contains information for efficient process- ing after the initial call and INBV must not be changed by the user. Distinct splines require distinct INBV parameters. WORK - work vector of length 3*K. DBVALU - value of the IDERIV-th derivative at X Error Conditions An improper input is a fatal error
DBVDER
SUBROUTINE DBVDER(X,Y,YP,G,IPAR)
DBVLUE
DOUBLE PRECISION FUNCTION DBVLUE(T,A,N,K,X,IDERIV) *** Double Precision version of BVALUE **** Calculates value at *X* of *IDERIV*-th deriv. of spline from B-repr.
DBVPOR
SUBROUTINE DBVPOR(Y,NROWY,NCOMP,XPTS,NXPTS,A,NROWA,ALPHA,NIC,B,
DCHFDV
SUBROUTINE DCHFDV(X1,X2,F1,F2,D1,D2,NE,XE,FE,DE,NEXT,IERR) **** Double Precision version of CHFDV **** DCHFDV: Cubic Hermite Function and Derivative Evaluator Evaluates the cubic polynomial determined by function values F1,F2 and derivatives D1,D2 on interval (X1,X2), together with its first derivative, at the points XE(J), J=1(1)NE. If only function values are required, use DCHFEV, instead.
Calling sequence: INTEGER NE, NEXT(2), IERR DOUBLE PRECISION X1, X2, F1, F2, D1, D2, XE(NE), FE(NE), DE(NE) CALL DCHFDV (X1,X2, F1,F2, D1,D2, NE, XE, FE, DE, NEXT, IERR) Parameters: X1,X2 -- (input) endpoints of interval of definition of cubic. (Error return if X1.EQ.X2 .) F1,F2 -- (input) values of function at X1 and X2, respectively. D1,D2 -- (input) values of derivative at X1 and X2, respectively. NE -- (input) number of evaluation points. (Error return if NE.LT.1 .) XE -- (input) real*8 array of points at which the functions are to be evaluated. If any of the XE are outside the interval [X1,X2], a warning error is returned in NEXT. FE -- (output) real*8 array of values of the cubic function defined by X1,X2, F1,F2, D1,D2 at the points XE. DE -- (output) real*8 array of values of the first derivative of the same function at the points XE. NEXT -- (output) integer array indicating number of extrapolation points: NEXT(1) = number of evaluation points to left of interval. NEXT(2) = number of evaluation points to right of interval. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable" errors: IERR = -1 if NE.LT.1 . IERR = -2 if X1.EQ.X2 . (Output arrays have not been changed in either case.)DCHFEV
SUBROUTINE DCHFEV(X1,X2,F1,F2,D1,D2,NE,XE,FE,NEXT,IERR) **** Double Precision version of CHFEV **** DCHFEV: Cubic Hermite Function EValuator Evaluates the cubic polynomial determined by function values F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points XE(J), J=1(1)NE.
Calling sequence: INTEGER NE, NEXT(2), IERR DOUBLE PRECISION X1, X2, F1, F2, D1, D2, XE(NE), FE(NE) CALL DCHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR) Parameters: X1,X2 -- (input) endpoints of interval of definition of cubic. (Error return if X1.EQ.X2 .) F1,F2 -- (input) values of function at X1 and X2, respectively. D1,D2 -- (input) values of derivative at X1 and X2, respectively. NE -- (input) number of evaluation points. (Error return if NE.LT.1 .) XE -- (input) real*8 array of points at which the function is to be evaluated. If any of the XE are outside the interval [X1,X2], a warning error is returned in NEXT. FE -- (output) real*8 array of values of the cubic function defined by X1,X2, F1,F2, D1,D2 at the points XE. NEXT -- (output) integer array indicating number of extrapolation points: NEXT(1) = number of evaluation points to left of interval. NEXT(2) = number of evaluation points to right of interval. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable" errors: IERR = -1 if NE.LT.1 . IERR = -2 if X1.EQ.X2 . (The FE-array has not been changed in either case.)DCHFIV
DOUBLE PRECISION FUNCTION DCHFIV(X1,X2,F1,F2,D1,D2,A,B,IERR) DCHFIV: Cubic Hermite Function Integral Evaluator. Called by DPCHIA to evaluate the integral of a single cubic (in Hermite form) over an arbitrary interval (A,B).
Calling sequence: INTEGER IERR DOUBLE PRECISION X1, X2, F1, F2, D1, D2, A, B DOUBLE PRECISION VALUE, DCHFIV VALUE = DCHFIV (X1, X2, F1, F2, D1, D2, A, B, IERR) Parameters: VALUE -- (output) VALUE of the requested integral. X1,X2 -- (input) endpoints if interval of definition of cubic. (Must be distinct. Error return if not.) F1,F2 -- (input) function values at the ends of the interval. D1,D2 -- (input) derivative values at the ends of the interval. A,B -- (input) endpoints of interval of integration. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable errors": IERR = -1 if X1.EQ.X2 . (VALUE has not been set in this case.)DCHFMC
INTEGER FUNCTION DCHFMC(D1,D2,DELTA) DCHFMC: Cubic Hermite Function Monotonicity Checker. Called by DPCHMC to determine the monotonicity properties of the cubic with boundary derivative values D1,D2 and chord slope DELTA.
Calling sequence: DOUBLE PRECISION D1, D2, DELTA INTEGER ISMON ISMON = DCHFMC (D1, D2, DELTA) Parameters: D1,D2 -- (input) derivative values at the ends of an interval. DELTA -- (input) data slope over that interval. ISMON -- (output) integer function value, indicating the monoto- nicity of the cubic segment: ISMON = -1 if function is strictly decreasing; ISMON = 0 if function is constant; ISMON = 1 if function is strictly increasing; ISMON = 2 if function is non-monotonic; ISMON = 3 if unable to determine. Fortran intrinsics used: DSIGN. Other routines used: D1MACH.DINTRV
SUBROUTINE DINTRV(XT,LXT,X,ILO,ILEFT,MFLAG) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June 1977, pp.441-472. Abstract **** a double precision routine **** DINTRV is the INTERV routine of the reference. DINTRV computes the largest integer ILEFT in 1 .LE. ILEFT .LE. LXT such that XT(ILEFT) .LE. X where XT(*) is a subdivision of the X interval. Precisely, X .LT. XT(1) 1 -1 if XT(I) .LE. X .LT. XT(I+1) then ILEFT=I , MFLAG=0 XT(LXT) .LE. X LXT 1, That is, when multiplicities are present in the break point to the left of X, the largest index is taken for ILEFT. Description of Arguments Input XT,X are double precision XT - XT is a knot or break point vector of length LXT LXT - length of the XT vector X - argument ILO - an initialization parameter which must be set to 1 the first time the spline array XT is processed by DINTRV. Output ILO - ILO contains information for efficient process- ing after the initial call and ILO must not be changed by the user. Distinct splines require distinct ILO parameters. ILEFT - largest integer satisfying XT(ILEFT) .LE. X MFLAG - signals when X lies out of bounds Error Conditions None
DINTYD
SUBROUTINE DINTYD(T,K,YH,NYH,DKY,IFLAG) DINTYD approximates the solution and derivatives at T by polynomial interpolation. Must be used in conjunction with the integrator package DDEBDF.
DINTYD computes interpolated values of the K-th derivative of the dependent variable vector Y, and stores it in DKY. This routine is called by DDEBDF with K = 0,1 and T = TOUT, but may also be called by the user for any K up to the current order. (see detailed instructions in LSODE usage documentation.)
The computed values in DKY are gotten by interpolation using the Nordsieck history array YH. This array corresponds uniquely to a vector-valued polynomial of degree NQCUR or less, and DKY is set to the K-th derivative of this polynomial at T. The formula for DKY is.. Q DKY(I) = Sum C(J,K) * (T - TN)**(J-K) * H**(-J) * YH(I,J+1) J=K where C(J,K) = J*(J-1)*...*(J-K+1), Q = NQCUR, TN = TCUR, H = HCUR. The quantities NQ = NQCUR, L = NQ+1, N = NEQ, TN, and H are communicated by common. The above sum is done in reverse order. IFLAG is returned negative if either K or T is out of bounds.
DJAIRY
SUBROUTINE DJAIRY(X,RX,C,AI,DAI)
DPCHFD
SUBROUTINE DPCHFD(N,X,F,D,INCFD,SKIP,NE,XE,FE,DE,IERR) **** Double Precision version of PCHFD **** DPCHFD: Piecewise Cubic Hermite Function and Derivative evaluator Evaluates the cubic Hermite function defined by N, X, F, D, to- gether with its first derivative, at the points XE(J), J=1(1)NE. If only function values are required, use DPCHFE, instead. To provide compatibility with DPCHIM and DPCHIC, includes an increment between successive values of the F- and D-arrays.
Calling sequence: PARAMETER (INCFD = ...) INTEGER N, NE, IERR DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE), DE(NE) LOGICAL SKIP CALL DPCHFD (N, X, F, D, INCFD, SKIP, NE, XE, FE, DE, IERR) Parameters: N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real*8 array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real*8 array of function values. F(1+(I-1)*INCFD) is the value corresponding to X(I). D -- (input) real*8 array of derivative values. D(1+(I-1)*INCFD) is the value corresponding to X(I). INCFD -- (input) increment between successive values in F and D. (Error return if INCFD.LT.1 .) SKIP -- (input/output) logical variable which should be set to .TRUE. if the user wishes to skip checks for validity of preceding parameters, or to .FALSE. otherwise. This will save time in case these checks have already been performed (say, in DPCHIM or DPCHIC). SKIP will be set to .TRUE. on normal return. NE -- (input) number of evaluation points. (Error return if NE.LT.1 .) XE -- (input) real*8 array of points at which the functions are to be evaluated. NOTES: 1. The evaluation will be most efficient if the elements of XE are increasing relative to X; that is, XE(J) .GE. X(I) implies XE(K) .GE. X(I), all K.GE.J . 2. If any of the XE are outside the interval [X(1),X(N)], values are extrapolated from the nearest extreme cubic, and a warning error is returned. FE -- (output) real*8 array of values of the cubic Hermite function defined by N, X, F, D at the points XE. DE -- (output) real*8 array of values of the first derivative of the same function at the points XE. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). Warning error: IERR.GT.0 means that extrapolation was performed at IERR points. "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. IERR = -4 if NE.LT.1 . (Output arrays have not been changed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated. IERR = -5 if an error has occurred in the lower-level routine DCHFDV. NB: this should never happen. Notify the author **IMMEDIATELY** if it does.DPCHFE
SUBROUTINE DPCHFE(N,X,F,D,INCFD,SKIP,NE,XE,FE,IERR) **** Double Precision version of PCHFE **** DPCHFE: Piecewise Cubic Hermite Function Evaluator Evaluates the cubic Hermite function defined by N, X, F, D at the points XE(J), J=1(1)NE. To provide compatibility with DPCHIM and DPCHIC, includes an increment between successive values of the F- and D-arrays.
Calling sequence: PARAMETER (INCFD = ...) INTEGER N, NE, IERR DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE) LOGICAL SKIP CALL DPCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR) Parameters: N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real*8 array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real*8 array of function values. F(1+(I-1)*INCFD) is the value corresponding to X(I). D -- (input) real*8 array of derivative values. D(1+(I-1)*INCFD) is the value corresponding to X(I). INCFD -- (input) increment between successive values in F and D. (Error return if INCFD.LT.1 .) SKIP -- (input/output) logical variable which should be set to .TRUE. if the user wishes to skip checks for validity of preceding parameters, or to .FALSE. otherwise. This will save time in case these checks have already been performed (say, in DPCHIM or DPCHIC). SKIP will be set to .TRUE. on normal return. NE -- (input) number of evaluation points. (Error return if NE.LT.1 .) XE -- (input) real*8 array of points at which the function is to be evaluated. NOTES: 1. The evaluation will be most efficient if the elements of XE are increasing relative to X; that is, XE(J) .GE. X(I) implies XE(K) .GE. X(I), all K.GE.J . 2. If any of the XE are outside the interval [X(1),X(N)], values are extrapolated from the nearest extreme cubic, and a warning error is returned. FE -- (output) real*8 array of values of the cubic Hermite function defined by N, X, F, D at the points XE. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). Warning error: IERR.GT.0 means that extrapolation was performed at IERR points. "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. IERR = -4 if NE.LT.1 . (The FE-array has not been changed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated.DPCHIA
DOUBLE PRECISION FUNCTION DPCHIA(N,X,F,D,INCFD,SKIP,A,B,IERR) **** Double Precision version of PCHIA **** DPCHIA: Piecewise Cubic Hermite Integrator, Arbitrary limits Evaluates the definite integral of the cubic Hermite function defined by N, X, F, D over the interval [A, B]. To provide compatibility with DPCHIM and DPCHIC, includes an increment between successive values of the F- and D-arrays.
Calling sequence: PARAMETER (INCFD = ...) INTEGER N, IERR DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N), A, B DOUBLE PRECISION VALUE, DPCHIA LOGICAL SKIP VALUE = DPCHIA (N, X, F, D, INCFD, SKIP, A, B, IERR) Parameters: VALUE -- (output) VALUE of the requested integral. N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real*8 array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real*8 array of function values. F(1+(I-1)*INCFD) is the value corresponding to X(I). D -- (input) real*8 array of derivative values. D(1+(I-1)*INCFD) is the value corresponding to X(I). INCFD -- (input) increment between successive values in F and D. (Error return if INCFD.LT.1 .) SKIP -- (input/output) logical variable which should be set to .TRUE. if the user wishes to skip checks for validity of preceding parameters, or to .FALSE. otherwise. This will save time in case these checks have already been performed (say, in DPCHIM or DPCHIC). SKIP will be set to .TRUE. on return with IERR.GE.0 . A,B -- (input) the limits of integration. NOTE: There is no requirement that [A,B] be contained in [X(1),X(N)]. However, the resulting integral value will be highly suspect, if not. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). Warning errors: IERR = 1 if A is outside the interval [X(1),X(N)]. IERR = 2 if B is outside the interval [X(1),X(N)]. IERR = 3 if both of the above are true. (Note that this means that either [A,B] contains data interval or the intervals do not intersect at all.) "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. (Value has not been computed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated.DPCHIC
SUBROUTINE DPCHIC(IC,VC,SWITCH,N,X,F,D,INCFD,WK,NWK,IERR) **** Double Precision version of PCHIC **** DPCHIC: Piecewise Cubic Hermite Interpolation Coefficients. Sets derivatives needed to determine a piecewise monotone piece- wise cubic interpolant to the data given in X and F satisfying the boundary conditions specified by IC and VC. The treatment of points where monotonicity switches direction is controlled by argument SWITCH. To facilitate two-dimensional applications, includes an increment between successive values of the F- and D-arrays. The resulting piecewise cubic Hermite function may be evaluated by DPCHFE or DPCHFD.
Calling sequence: PARAMETER (INCFD = ...) INTEGER IC(2), N, NWK, IERR DOUBLE PRECISION VC(2), SWITCH, X(N), F(INCFD,N), D(INCFD,N), WK(NWK) CALL DPCHIC (IC, VC, SWITCH, N, X, F, D, INCFD, WK, NWK, IERR) Parameters: IC -- (input) integer array of length 2 specifying desired boundary conditions: IC(1) = IBEG, desired condition at beginning of data. IC(2) = IEND, desired condition at end of data. IBEG = 0 for the default boundary condition (the same as used by DPCHIM). If IBEG.NE.0, then its sign indicates whether the boundary derivative is to be adjusted, if necessary, to be compatible with monotonicity: IBEG.GT.0 if no adjustment is to be performed. IBEG.LT.0 if the derivative is to be adjusted for monotonicity. Allowable values for the magnitude of IBEG are: IBEG = 1 if first derivative at X(1) is given in VC(1). IBEG = 2 if second derivative at X(1) is given in VC(1). IBEG = 3 to use the 3-point difference formula for D(1). (Reverts to the default b.c. if N.LT.3 .) IBEG = 4 to use the 4-point difference formula for D(1). (Reverts to the default b.c. if N.LT.4 .) IBEG = 5 to set D(1) so that the second derivative is con- tinuous at X(2). (Reverts to the default b.c. if N.LT.4.) This option is somewhat analogous to the "not a knot" boundary condition provided by DPCHSP. NOTES (IBEG): 1. An error return is taken if DABS(IBEG).GT.5 . 2. Only in case IBEG.LE.0 is it guaranteed that the interpolant will be monotonic in the first interval. If the returned value of D(1) lies between zero and 3*SLOPE(1), the interpolant will be monotonic. This is **NOT** checked if IBEG.GT.0 . 3. If IBEG.LT.0 and D(1) had to be changed to achieve mono- tonicity, a warning error is returned. IEND may take on the same values as IBEG, but applied to derivative at X(N). In case IEND = 1 or 2, the value is given in VC(2). NOTES (IEND): 1. An error return is taken if DABS(IEND).GT.5 . 2. Only in case IEND.LE.0 is it guaranteed that the interpolant will be monotonic in the last interval. If the returned value of D(1+(N-1)*INCFD) lies between zero and 3*SLOPE(N-1), the interpolant will be monotonic. This is **NOT** checked if IEND.GT.0 . 3. If IEND.LT.0 and D(1+(N-1)*INCFD) had to be changed to achieve monotonicity, a warning error is returned. VC -- (input) real*8 array of length 2 specifying desired boundary values, as indicated above. VC(1) need be set only if IC(1) = 1 or 2 . VC(2) need be set only if IC(2) = 1 or 2 . SWITCH -- (input) indicates desired treatment of points where direction of monotonicity switches: Set SWITCH to zero if interpolant is required to be mono- tonic in each interval, regardless of monotonicity of data. NOTES: 1. This will cause D to be set to zero at all switch points, thus forcing extrema there. 2. The result of using this option with the default boun- dary conditions will be identical to using DPCHIM, but will generally cost more compute time. This option is provided only to facilitate comparison of different switch and/or boundary conditions. Set SWITCH nonzero to use a formula based on the 3-point difference formula in the vicinity of switch points. If SWITCH is positive, the interpolant on each interval containing an extremum is controlled to not deviate from the data by more than SWITCH*DFLOC, where DFLOC is the maximum of the change of F on this interval and its two immediate neighbors. If SWITCH is negative, no such control is to be imposed. N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real*8 array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real*8 array of dependent variable values to be interpolated. F(1+(I-1)*INCFD) is value corresponding to X(I). D -- (output) real*8 array of derivative values at the data points. These values will determine a monotone cubic Hermite function on each subinterval on which the data are monotonic, except possibly adjacent to switches in monotonicity. The value corresponding to X(I) is stored in D(1+(I-1)*INCFD), I=1(1)N. No other entries in D are changed. INCFD -- (input) increment between successive values in F and D. This argument is provided primarily for 2-D applications. (Error return if INCFD.LT.1 .) WK -- (scratch) real*8 array of working storage. The user may wish to know that the returned values are: WK(I) = H(I) = X(I+1) - X(I) ; WK(N-1+I) = SLOPE(I) = (F(1,I+1) - F(1,I)) / H(I) for I = 1(1)N-1. NWK -- (input) length of work array. (Error return if NWK.LT.2*(N-1) .) IERR -- (output) error flag. Normal return: IERR = 0 (no errors). Warning errors: IERR = 1 if IBEG.LT.0 and D(1) had to be adjusted for monotonicity. IERR = 2 if IEND.LT.0 and D(1+(N-1)*INCFD) had to be adjusted for monotonicity. IERR = 3 if both of the above are true. "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. IERR = -4 if DABS(IBEG).GT.5 . IERR = -5 if DABS(IEND).GT.5 . IERR = -6 if both of the above are true. IERR = -7 if NWK.LT.2*(N-1) . (The D-array has not been changed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated.DPCHID
DOUBLE PRECISION FUNCTION DPCHID(N,X,F,D,INCFD,SKIP,IA,IB,IERR) **** Double Precision version of PCHID **** DPCHID: Piecewise Cubic Hermite Integrator, Data limits Evaluates the definite integral of the cubic Hermite function defined by N, X, F, D over the interval [X(IA), X(IB)]. To provide compatibility with DPCHIM and DPCHIC, includes an increment between successive values of the F- and D-arrays.
Calling sequence: PARAMETER (INCFD = ...) INTEGER N, IA, IB, IERR DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N) LOGICAL SKIP VALUE = DPCHID (N, X, F, D, INCFD, SKIP, IA, IB, IERR) Parameters: VALUE -- (output) VALUE of the requested integral. N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real*8 array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real*8 array of function values. F(1+(I-1)*INCFD) is the value corresponding to X(I). D -- (input) real*8 array of derivative values. D(1+(I-1)*INCFD) is the value corresponding to X(I). INCFD -- (input) increment between successive values in F and D. (Error return if INCFD.LT.1 .) SKIP -- (input/output) logical variable which should be set to .TRUE. if the user wishes to skip checks for validity of preceding parameters, or to .FALSE. otherwise. This will save time in case these checks have already been performed (say, in DPCHIM or DPCHIC). SKIP will be set to .TRUE. on return with IERR = 0 or -4. IA,IB -- (input) indices in X-array for the limits of integration. both must be in the range [1,N]. (Error return if not.) No restrictions on their relative values. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. IERR = -4 if IA or IB is out of range. (Value has not been computed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated.DPCHIM
SUBROUTINE DPCHIM(N,X,F,D,INCFD,IERR) **** Double Precision version of PCHIM **** DPCHIM: Piecewise Cubic Hermite Interpolation to Monotone data. Sets derivatives needed to determine a monotone piecewise cubic Hermite interpolant to the data given in X and F. Default boundary conditions are provided which are compatible with monotonicity. (See DPCHIC if user control of boundary con- ditions is desired.) If the data are only piecewise monotonic, the interpolant will have an extremum at each point where monotonicity switches direc- tion. (See DPCHIC if user control is desired in such cases.) To facilitate two-dimensional applications, includes an increment between successive values of the F- and D-arrays. The resulting piecewise cubic Hermite function may be evaluated by DPCHFE or DPCHFD.
Calling sequence: PARAMETER (INCFD = ...) INTEGER N, IERR DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N) CALL DPCHIM (N, X, F, D, INCFD, IERR) Parameters: N -- (input) number of data points. (Error return if N.LT.2 .) If N=2, simply does linear interpolation. X -- (input) real*8 array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real*8 array of dependent variable values to be interpolated. F(1+(I-1)*INCFD) is value corresponding to X(I). DPCHIM is designed for monotonic data, but it will work for any F-array. It will force extrema at points where monotonicity switches direction. If some other treatment of switch points is desired, DPCHIC should be used instead. ----- D -- (output) real*8 array of derivative values at the data points. If the data are monotonic, these values will determine a monotone cubic Hermite function. The value corresponding to X(I) is stored in D(1+(I-1)*INCFD), I=1(1)N. No other entries in D are changed. INCFD -- (input) increment between successive values in F and D. This argument is provided primarily for 2-D applications. (Error return if INCFD.LT.1 .) IERR -- (output) error flag. Normal return: IERR = 0 (no errors). Warning error: IERR.GT.0 means that IERR switches in the direction of monotonicity were detected. "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. (The D-array has not been changed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated.DPCHMC
SUBROUTINE DPCHMC(N,X,F,D,INCFD,SKIP,ISMON,IERR) **** Double Precision version of PCHMC **** DPCHMC: Piecewise Cubic Hermite Monotonicity Checker. Checks the cubic Hermite function defined by N, X, F, D for monotonicity. To provide compatibility with DPCHIM and DPCHIC, includes an increment between successive values of the F- and D-arrays.
Calling sequence: PARAMETER (INCFD = ...) INTEGER N, ISMON(N), IERR DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N) LOGICAL SKIP CALL DPCHMC (N, X, F, D, INCFD, SKIP, ISMON, IERR) Parameters: N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real*8 array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real*8 array of function values. F(1+(I-1)*INCFD) is the value corresponding to X(I). D -- (input) real*8 array of derivative values. D(1+(I-1)*INCFD) is the value corresponding to X(I). INCFD -- (input) increment between successive values in F and D. (Error return if INCFD.LT.1 .) SKIP -- (input/output) logical variable which should be set to .TRUE. if the user wishes to skip checks for validity of preceding parameters, or to .FALSE. otherwise. This will save time in case these checks have already been performed. SKIP will be set to .TRUE. on normal return. ISMON -- (output) integer array indicating on which intervals the PCH function defined by N, X, F, D is monotonic. For data interval [X(I),X(I+1)], ISMON(I) = -1 if function is strictly decreasing; ISMON(I) = 0 if function is constant; ISMON(I) = 1 if function is strictly increasing; ISMON(I) = 2 if function is non-monotonic; ISMON(I) = 3 if unable to determine. (This means that the D-values are near the boundary of the monotonicity region. A small increase pro- duces non-monotonicity; decrease, strict monotonicity.) The above applies to I=1(1)N-1. ISMON(N) indicates whether the entire function is monotonic on [X(1),X(N)]. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. (The ISMON-array has not been changed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated.DPCHNG
SUBROUTINE DPCHNG(II,XVAL,IPLACE,SX,IX,IRCX) THIS SUBROUTINE IS A MODIFICATION OF THE SUBROUTINE LCHNGS, SANDIA LABS. REPT. SAND78-0785. DPCHNG LIMITS THE TYPE OF STORAGE TO A SEQUENTIAL SCHEME. MODIFICATIONS BY K.L. HIEBERT AND R.J. HANSON REVISED 811130-1000 REVISED YYMMDD-HHMM SPARSE MATRIX ELEMENT ALTERATION SUBROUTINE. AUTHORS JOHN A. WISNIEWSKI AND R.J. HANSON SUBROUTINE DPCHNG() CHANGES ELEMENT II IN VECTOR +/- IRCX TO THE VALUE XVAL. II THE ABSOLUTE VALUE OF THIS INTEGER IS THE SUBSCRIPT FOR THE ELEMENT TO BE CHANGED. XVAL NEW VALUE OF THE MATRIX ELEMENT BEING CHANGED. IPLACE POINTER INFORMATION WHICH IS MAINTAINED BY THE PACKAGE. SX(*),IX(*) THE WORK ARRAYS WHICH ARE USED TO STORE THE SPARSE MATRIX. THESE ARRAYS ARE AUTOMATICALLY MAINTAINED BY THE PACKAGE FOR THE USER. IRCX POINTS TO THE VECTOR OF THE MATRIX BEING UPDATED. A NEGATIVE VALUE OF IRCX INDICATES THAT ROW -IRCX IS BEING UPDATED. A POSITIVE VALUE OF IRCX INDICATES THAT COLUMN IRCX IS BEING UPDATED. A ZERO VALUE OF IRCX IS AN ERROR. SINCE DATA ITEMS ARE KEPT SORTED IN THE SEQUENTIAL DATA STRUCTURE, CHANGING A MATRIX ELEMENT CAN REQUIRE THE MOVEMENT OF ALL THE DATA ITEMS IN THE MATRIX. FOR THIS REASON, IT IS SUGGESTED THAT DATA ITEMS BE ADDED A COL. AT A TIME, IN ASCENDING COL. SEQUENCE. FURTHERMORE, SINCE DELETING ITEMS FROM THE DATA STRUCTURE MAY ALSO REQUIRE MOVING LARGE AMOUNTS OF DATA, ZERO ELEMENTS ARE EXPLICITLY STORED IN THE MATRIX.
DPCHSP
SUBROUTINE DPCHSP(IC,VC,N,X,F,D,INCFD,WK,NWK,IERR) **** Double Precision version of PCHSP **** DPCHSP: Piecewise Cubic Hermite Spline Computes the Hermite representation of the cubic spline inter- polant to the data given in X and F satisfying the boundary conditions specified by IC and VC. To facilitate two-dimensional applications, includes an increment between successive values of the F- and D-arrays. The resulting piecewise cubic Hermite function may be evaluated by DPCHFE or DPCHFD. NOTE: This is a modified version of C. de Boor'S cubic spline routine CUBSPL.
Calling sequence: PARAMETER (INCFD = ...) INTEGER IC(2), N, NWK, IERR DOUBLE PRECISION VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK) CALL DPCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR) Parameters: IC -- (input) integer array of length 2 specifying desired boundary conditions: IC(1) = IBEG, desired condition at beginning of data. IC(2) = IEND, desired condition at end of data. IBEG = 0 to set D(1) so that the third derivative is con- tinuous at X(2). This is the "not a knot" condition provided by de Boor'S cubic spline routine CUBSPL. < This is the default boundary condition. > IBEG = 1 if first derivative at X(1) is given in VC(1). IBEG = 2 if second derivative at X(1) is given in VC(1). IBEG = 3 to use the 3-point difference formula for D(1). (Reverts to the default b.c. if N.LT.3 .) IBEG = 4 to use the 4-point difference formula for D(1). (Reverts to the default b.c. if N.LT.4 .) NOTES: 1. An error return is taken if IBEG is out of range. 2. For the "natural" boundary condition, use IBEG=2 and VC(1)=0. IEND may take on the same values as IBEG, but applied to derivative at X(N). In case IEND = 1 or 2, the value is given in VC(2). NOTES: 1. An error return is taken if IEND is out of range. 2. For the "natural" boundary condition, use IEND=2 and VC(2)=0. VC -- (input) real*8 array of length 2 specifying desired boundary values, as indicated above. VC(1) need be set only if IC(1) = 1 or 2 . VC(2) need be set only if IC(2) = 1 or 2 . N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real*8 array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real*8 array of dependent variable values to be interpolated. F(1+(I-1)*INCFD) is value corresponding to X(I). D -- (output) real*8 array of derivative values at the data points. These values will determine the cubic spline interpolant with the requested boundary conditions. The value corresponding to X(I) is stored in D(1+(I-1)*INCFD), I=1(1)N. No other entries in D are changed. INCFD -- (input) increment between successive values in F and D. This argument is provided primarily for 2-D applications. (Error return if INCFD.LT.1 .) WK -- (scratch) real*8 array of working storage. NWK -- (input) length of work array. (Error return if NWK.LT.2*N .) IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. IERR = -4 if IBEG.LT.0 or IBEG.GT.4 . IERR = -5 if IEND.LT.0 of IEND.GT.4 . IERR = -6 if both of the above are true. IERR = -7 if NWK is too small. NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated. (The D-array has not been changed in any of these cases.) IERR = -8 in case of trouble solving the linear system for the interior derivative values. (The D-array may have been changed in this case.) ( Do **NOT** use it! )DPCHST
DOUBLE PRECISION FUNCTION DPCHST(ARG1,ARG2) DPCHST: DPCHIP Sign-Testing Routine. Returns: -1. if ARG1 and ARG2 are of opposite sign. 0. if either argument is zero. +1. if ARG1 and ARG2 are of the same sign. The object is to do this without multiplying ARG1*ARG2, to avoid possible over/underflow problems. Fortran intrinsics used: SIGN.
Programmed by: Fred N. Fritsch, FTS 532-4275, (415) 422-4275, Mathematics and Statistics Division, Lawrence Livermore National Laboratory. Change record: 82-08-05 Converted to SLATEC library version.
Programming notes: To produce a single precision version, simply: a. Change DPCHST to PCHST wherever it occurs, b. Change all references to the Fortran intrinsics to their single presision equivalents, c. Change the double precision declarations to real, and d. Change the constants ZERO and ONE to single precision.DPCHSW
SUBROUTINE DPCHSW(DFMAX,IEXTRM,D1,D2,H,SLOPE,IERR) DPCHSW: DPCHCS Switch Excursion Limiter. Called by DPCHCS to adjust D1 and D2 if necessary to insure that the extremum on this interval is not further than DFMAX from the extreme data value.
Calling sequence: INTEGER IEXTRM, IERR DOUBLE PRECISION DFMAX, D1, D2, H, SLOPE CALL DPCHSW (DFMAX, IEXTRM, D1, D2, H, SLOPE, IERR) Parameters: DFMAX -- (input) maximum allowed difference between F(IEXTRM) and the cubic determined by derivative values D1,D2. (assumes DFMAX.GT.0.) IEXTRM -- (input) index of the extreme data value. (assumes IEXTRM = 1 or 2 . Any value .NE.1 is treated as 2.) D1,D2 -- (input) derivative values at the ends of the interval. (Assumes D1*D2 .LE. 0.) (output) may be modified if necessary to meet the restriction imposed by DFMAX. H -- (input) interval length. (Assumes H.GT.0.) SLOPE -- (input) data SLOPE on the interval. IERR -- (output) error flag. should be zero. If IERR=-1, assumption on D1 and D2 is not satisfied. If IERR=-2, quadratic equation locating extremum has negative descriminant (should never occur). ------- WARNING: This routine does no validity-checking of arguments. ------- Fortran intrinsics used: DABS, SIGN.DPLINT
SUBROUTINE DPLINT(N,X,Y,C) *** Double Precision Version of POLINT *** Written by Robert E. Huddleston, Sandia Laboratories, Livermore Abstract Subroutine DPLINT is designed to produce the polynomial which interpolates the data (X(I),Y(I)), I=1,...,N. DPLINT sets up information in the array C which can be used by subroutine DPOLVL to evaluate the polynomial and its derivatives and by subroutine DPOLCF to produce the coefficients. Formal Parameters *** All TYPE REAL variables are DOUBLE PRECISION *** N - the number of data points (N .GE. 1) X - the array of abscissas (all of which must be distinct) Y - the array of ordinates C - an array of information used by subroutines ******* Dimensioning Information ******* Arrays X,Y, and C must be dimensioned at least N in the calling program.
DPLPCE
SUBROUTINE DPLPCE(MRELAS,NVARS,LMX,LBM,ITLP,ITBRC,IBASIS,IMAT, THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. /REAL (12 BLANKS)/DOUBLE PRECISION/,/ABS(/DABS(/,/AMAX1/ DMAX1/,/SASUM/DASUM/,/DCOPY/,DCOPY/. REVISED 811219-1630 REVISED YYMMDD-HHMM THIS SUBPROGRAM IS FROM THE DSPLP( ) PACKAGE. IT CALCULATES THE APPROXIMATE ERROR IN THE PRIMAL AND DUAL SYSTEMS. IT IS THE MAIN PART OF THE PROCEDURE (COMPUTE ERROR IN DUAL AND PRIMAL SYSTEMS).
DPLPDM
SUBROUTINE DPLPDM(MRELAS,NVARS,LMX,LBM,NREDC,INFO,IOPT,IBASIS, THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. /REAL (12 BLANKS)/DOUBLE PRECISION/. THIS SUBPROGRAM IS FROM THE DSPLP( ) PACKAGE. IT PERFORMS THE TASK OF DEFINING THE ENTRIES OF THE BASIS MATRIX AND DECOMPOSING IT USING THE LA05 PACKAGE. IT IS THE MAIN PART OF THE PROCEDURE (DECOMPOSE BASIS MATRIX). REVISED 811123-1500 REVISED YYMMDD-HHMM
DPLPFE
SUBROUTINE DPLPFE(MRELAS,NVARS,LMX,LBM,IENTER,IBASIS,IMAT,IBRC, THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. /REAL (12 BLANKS)/DOUBLE PRECISION/,/ABS(/DABS(/,/SASUM/DASUM/, /SCOPY/DCOPY/. THIS SUBPROGRAM IS PART OF THE DSPLP( ) PACKAGE. IT IMPLEMENTS THE PROCEDURE (FIND VARIABLE TO ENTER BASIS AND GET SEARCH DIRECTION). REVISED 811130-1100 REVISED YYMMDD-HHMM
DPLPFL
SUBROUTINE DPLPFL(MRELAS,NVARS,IENTER,ILEAVE,IBASIS,IND,IBB, THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. /REAL (12 BLANKS)/DOUBLE PRECISION/,/ABS(/DABS(/. THIS SUBPROGRAM IS PART OF THE DSPLP( ) PACKAGE. IT IMPLEMENTS THE PROCEDURE (CHOOSE VARIABLE TO LEAVE BASIS). REVISED 811130-1045 REVISED YYMMDD-HHMM
DPLPMN
SUBROUTINE DPLPMN(DUSRMT,MRELAS,NVARS,COSTS,PRGOPT,DATTRV,BL,BU, THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. REAL (12 BLANKS)/DOUBLE PRECISION/,/SASUM/DASUM/,/SDOT/DDOT/, /SCOPY/DCOPY/,/ABS(/DABS(/,/AMAX1/DMAX1/,/AMIN1/DMIN1/, /SVOUT/DVOUT/. MARVEL OPTION(S).. OUTPUT=YES/NO TO ELIMINATE PRINTED OUTPUT. THIS DOES NOT APPLY TO THE CALLS TO THE ERROR PROCESSOR. MAIN SUBROUTINE FOR DSPLP(*) PACKAGE.
DPLPMU
SUBROUTINE DPLPMU(MRELAS,NVARS,LMX,LBM,NREDC,INFO,IENTER,ILEAVE, THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. /REAL (12 BLANKS)/DOUBLE PRECISION/,/ABS(/DABS(/,/AMAX1/DMAX1/, /SASUM/DASUM/,/SCOPY/DCOPY/,/SDOT/DDOT/, /.E0/.D0/ THIS SUBPROGRAM IS FROM THE DSPLP( ) PACKAGE. IT PERFORMS THE TASKS OF UPDATING THE PRIMAL SOLUTION, EDGE WEIGHTS, REDUCED COSTS, AND MATRIX DECOMPOSITION. IT IS THE MAIN PART OF THE PROCEDURE (MAKE MOVE AND UPDATE). REVISED 821122-1100 REVISED YYMMDD
DPLPUP
SUBROUTINE DPLPUP(DUSRMT,MRELAS,NVARS,PRGOPT,DATTRV,BL,BU,IND, THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. /REAL (12 BLANKS)/DOUBLE PRECISION/,/ABS(/DABS(/. REVISED 810613-1130 REVISED YYMMDD-HHMM THIS SUBROUTINE COLLECTS INFORMATION ABOUT THE BOUNDS AND MATRIX FROM THE USER. IT IS PART OF THE DSPLP( ) PACKAGE.
DPNNZR
SUBROUTINE DPNNZR(I,XVAL,IPLACE,SX,IX,IRCX) THIS SUBROUTINE IS A MODIFICATION OF THE SUBROUTINE LNNZRS, SANDIA LABS. REPT. SAND78-0785. DPNNZR LIMITS THE TYPE OF STORAGE TO A SEQUENTIAL SCHEME. MODIFICATIONS BY K.L. HIEBERT AND R.J. HANSON REVISED 811130-1000 REVISED YYMMDD-HHMM SPARSE MATRIX NON ZERO RETRIEVAL SUBROUTINE. AUTHORS JOHN A. WISNIEWSKI AND R. J. HANSON SUBROUTINE DPNNZR() GETS THE NEXT NONZERO VALUE IN ROW OR COLUMN +/- IRCX WITH AN INDEX GREATER THAN THE VALUE OF I. I ABSOLUTE VALUE OF THIS SUBSCRIPT IS TO BE EXCEEDED IN THE SEARCH FOR THE NEXT NONZERO VALUE. A NEGATIVE OR ZERO VALUE OF I CAUSES THE SEARCH TO START AT THE BEGINNING OF THE VECTOR. A POSITIVE VALUE OF I CAUSES THE SEARCH TO CONTINUE FROM THE LAST PLACE ACCESSED. ON OUTPUT, THE ARGUMENT I CONTAINS THE VALUE OF THE SUBSCRIPT FOUND. AN OUTPUT VALUE OF I EQUAL TO ZERO INDICATES THAT ALL COMPONENTS WITH AN INDEX GREATER THAN THE INPUT VALUE OF I ARE ZERO. XVAL VALUE OF THE NONZERO ELEMENT FOUND. ON OUTPUT, XVAL=0. WHENEVER I=0. IPLACE POINTER INFORMATION WHICH IS MAINTAINED BY THE PACKAGE. SX(*),IX(*) THE WORK ARRAYS WHICH ARE USED TO STORE THE SPARSE MATRIX. THESE ARRAY CONTENTS ARE AUTOMATICALLY MAINTAINED BY THE PACKAGE FOR THE USER. IRCX POINTS TO THE VECTOR OF THE MATRIX BEING SCANNED. A NEGATIVE VALUE OF IRCX INDICATES THAT ROW -IRCX IS TO BE SCANNED. A POSITIVE VALUE OF IRCX INDICATES THAT COLUMN IRCX IS TO BE SCANNED. A ZERO VALUE OF IRCX IS AN ERROR.
DPOLVL
SUBROUTINE DPOLVL(NDER,XX,YFIT,YP,N,X,C,WORK,IERR) **** Double Precision version of POLYVL **** Written by Robert E. Huddleston, Sandia Laboratories, Livermore Abstract - Subroutine DPOLVL calculates the value of the polynomial and its first NDER derivatives where the polynomial was produced by a previous call to HRMITE or DPLINT. The variable N and the arrays X and C must not be altered between the call to HRMITE or DPLINT and the call to DPOLVL. ****** Dimensioning Information ******* YP must be dimensioned by at least NDER X must be dimensioned by at least N (see the abstract ) C must be dimensioned by at least N (see the abstract ) WORK must be dimensioned by at least 2*N if NDER is .GT. 0. *** Note *** If NDER=0, neither YP nor WORK need to be dimensioned variables. If NDER=1, YP does not need to be a dimensioned variable. ***** Input parameters *** All TYPE REAL variables are DOUBLE PRECISION *** NDER - the number of derivatives to be evaluated XX - the argument at which the polynomial and its derivatives are to be evaluated. N - ***** * N, X, and C must not be altered between the call X - * to HRMITE (or the call to DPLINT) and the call * to DPOLVL. C - ***** ***** Output Parameters *** All TYPE REAL variables are DOUBLE PRECISION *** YFIT - the value of the polynomial at XX YP - the derivatives of the polynomial at XX. The derivative of order J at XX is stored in YP(J) , J = 1,...,NDER. IERR - Output error flag with the following possible values. = 1 indicates normal execution ***** Storage Parameters WORK = this is an array to provide internal working storage for DPOLVL. It must be dimensioned by at least 2*N if NDER is .GT. 0. If NDER=0, WORK does not need to be a dimensioned variable.
DPOPT
SUBROUTINE DPOPT(PRGOPT,MRELAS,NVARS,INFO,CSC,IBASIS,ROPT,INTOPT, THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES. USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/. /REAL (12 BLANKS)/DOUBLE PRECISION/,/ABS(/DABS(/,/AMAX1/ DMAX1/,SQRT/DSQRT/,/R1MACH/D1MACH/,/AMIN1/DMIN1/,/E0/D0/ REVISED 821122-1045 REVISED YYMMDD-HHMM THIS SUBROUTINE PROCESSES THE OPTION VECTOR, PRGOPT(*), AND VALIDATES ANY MODIFIED DATA.
DPPVAL
DOUBLE PRECISION FUNCTION DPPVAL(LDC,C,XI,LXI,K,IDERIV,X,INPPV) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract **** a double precision routine **** DPPVAL is the PPVALU function of the reference. DPPVAL calculates (at X) the value of the IDERIV-th derivative of the B-spline from the PP-representation (C,XI,LXI,K). The Taylor expansion about XI(J) for X in the interval XI(J) .LE. X .LT. XI(J+1) is evaluated, J=1,LXI. Right limiting values at X=XI(J) are obtained. DPPVAL will extrapolate beyond XI(1) and XI(LXI+1). To obtain left limiting values (left derivatives) at XI(J) replace LXI by J-1 and set X=XI(J),J=2,LXI+1. DPPVAL calls DINTRV Description of Arguments Input C,XI,X are double precision LDC - leading dimension of C matrix, LDC .GE. K C - matrix of dimension at least (K,LXI) containing right derivatives at break points XI(*). XI - break point vector of length LXI+1 LXI - number of polynomial pieces K - order of B-spline, K .GE. 1 IDERIV - order of the derivative, 0 .LE. IDERIV .LE. K-1 IDERIV=0 gives the B-spline value X - argument, XI(1) .LE. X .LE. XI(LXI+1) INPPV - an initialization parameter which must be set to 1 the first time DPPVAL is called. Output DPPVAL is double precision INPPV - INPPV contains information for efficient process- ing after the initial call and INPPV must not be changed by the user. Distinct splines require distinct INPPV parameters. DPPVAL - value of the IDERIV-th derivative at X Error Conditions Improper input is a fatal error
DPRVEC
DOUBLE PRECISION FUNCTION DPRVEC(M,U,V)
DPRWPG
SUBROUTINE DPRWPG(KEY,IPAGE,LPG,SX,IX) THIS SUBROUTINE IS A MODIFICATION OF THE SUBROUTINE LRWPGE, SANDIA LABS. REPT. SAND78-0785. DPRWPG LIMITS THE TYPE OF STORAGE TO A SEQUENTIAL SCHEME. MODIFICATIONS BY K.L. HIEBERT AND R.J. HANSON REVISED 811130-1000 REVISED YYMMDD-HHMM VIRTUAL MEMORY PAGE READ/WRITE SUBROUTINE. AUTHORS J. A. WISNIEWSKI AND R.J. HANSON DEPENDING ON THE VALUE OF KEY, SUBROUTINE DPRWPG() PERFORMS A PAGE READ OR WRITE OF PAGE IPAGE. THE PAGE HAS LENGTH LPG. KEY IS A FLAG INDICATING WHETHER A PAGE READ OR WRITE IS TO BE PERFORMED. IF KEY = 1 DATA IS READ. IF KEY = 2 DATA IS WRITTEN. IPAGE IS THE PAGE NUMBER OF THE MATRIX TO BE ACCESSED. LPG IS THE LENGTH OF THE PAGE OF THE MATRIX TO BE ACCESSED. SX(*),IX(*) IS THE MATRIX TO BE ACCESSED. EXTERNAL REFERENCES XERROR,DPRWVR
DPRWVR
SUBROUTINE DPRWVR(KEY,IPAGE,LPG,SX,IX) THIS SUBROUTINE IS A MODIFICATION OF THE SUBROUTINE LRWVIR. DPRWVR LIMITS THE TYPE OF STORAGE TO A SEQUENTIAL SPARSE MATRIX STORAGE SCHEME. THE PAGE STORAGE IS ON RANDOM ACCESS DISK. DPRWVR() IS PART OF THE SPARSE LP PACKAGE, DSPLP( ). MODIFICATIONS BY K.L. HIEBERT AND R.J. HANSON AUTHORS J. A. WISNIEWSKI AND R.J. HANSON KEY IS A FLAG WHICH INDICATES WHETHER A READ OR WRITE OPERATION IS TO BE PERFORMED. A VALUE OF KEY=1 INDICATES A READ. A VALUE OF KEY=2 INDICATES A WRITE. IPAGE IS THE PAGE OF MATRIX MN WE ARE ACCESSING. LPG IS THE LENGTH OF THE PAGE. SX(*),IX(*) IS THE MATRIX DATA.
INTRV
SUBROUTINE INTRV(XT,LXT,X,ILO,ILEFT,MFLAG) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June 1977, pp. 441-472. Abstract INTRV is the INTERV routine of the reference. INTRV computes the largest integer ILEFT in 1 .LE. ILEFT .LE. LXT such that XT(ILEFT) .LE. X where XT(*) is a subdivision of the X interval. Precisely, X .LT. XT(1) 1 -1 if XT(I) .LE. X .LT. XT(I+1) then ILEFT=I , MFLAG=0 XT(LXT) .LE. X LXT 1, That is, when multiplicities are present in the break point to the left of X, the largest index is taken for ILEFT. Description of Arguments Input XT - XT is a knot or break point vector of length LXT LXT - length of the XT vector X - argument ILO - an initialization parameter which must be set to 1 the first time the spline array XT is processed by INTRV. Output ILO - ILO contains information for efficient process- ing after the initial call, and ILO must not be changed by the user. Distinct splines require distinct ILO parameters. ILEFT - largest integer satisfying XT(ILEFT) .LE. X MFLAG - signals when X lies out of bounds Error Conditions None
INTYD
SUBROUTINE INTYD(T,K,YH,NYH,DKY,IFLAG)
PCHDOC
SUBROUTINE PCHDOC PCHIP: Piecewise Cubic Hermite Interpolation Package This document contains the specifications for PCHIP, a new Fortran package for piecewise cubic Hermite interpolation of data. It features software to produce a monotone and "visually pleasing" interpolant to monotone data. As is demonstrated in Reference 1, such an interpolant may be more reasonable than a cubic spline if the data contains both "steep" and "flat" sections. Interpola- tion of cumulative probability distribution functions is another application. (See References 1-3 for examples.) All piecewise cubic functions in PCHIP are represented in cubic Hermite form; that is, F(X) is determined by its values F(I) and derivatives D(I) at the breakpoints X(I), I=1(1)N. The double precision equivalents of the PCHIP routines are obtained from the single precision names by prefixing the single precision names with a D. For example, the double precision equivalent of PCHIM is DPCHIM. The contents of the package are as follows: 1. Determine Derivative Values. PCHIM -- Piecewise Cubic Hermite Interpolation to Monotone data. Used if the data are monotonic or if the user wants to guarantee that the interpolant stays within the limits of the data. (See Reference 2.) PCHIC -- Piecewise Cubic Hermite Interpolation Coefficients. Used if neither of the above conditions holds, or if the user wishes control over boundary derivatives. Will generally reproduce monotonicity on subintervals over which the data are monotonic. PCHSP -- Piecewise Cubic Hermite Spline. Produces a cubic spline interpolator in cubic Hermite form. Provided primarily for easy comparison of the spline with other piecewise cubic interpolants. (A modified version of de Boor'S CUBSPL, Reference 4.) 2. Evaluate, Differentiate, or Integrate Resulting PCH Function. NOTE: If derivative values are available from some other source, these routines can be used without calling any of the previous routines. CHFEV -- Cubic Hermite Function EValuator. Evaluates a single cubic Hermite function at an array of points. Used when the interval is known, as in graphing applications. Called by PCHFE. PCHFE -- Piecewise Cubic Hermite Function Evaluator. Used when the interval is unknown or the evaluation array spans more than one data interval. CHFDV -- Cubic Hermite Function and Derivative Evaluator. Evaluates a single cubic Hermite function and its first derivative at an array of points. Used when the interval is known, as in graphing applications. Called by PCHFD. PCHFD -- Piecewise Cubic Hermite Function and Derivative Evaluator. Used when the interval is unknown or the evaluation array spans more than one data interval. PCHID -- Piecewise Cubic Hermite Integrator, Data Limits. Computes the definite integral of a piecewise cubic Hermite function when the integration limits are data points. PCHIA -- Piecewise Cubic Hermite Integrator, Arbitrary Limits. Computes the definite integral of a piecewise cubic Hermite function over an arbitrary finite interval. 3. Check for monotonicity. PCHMC -- Piecewise Cubic Hermite Monotonicity Checker. 4. Internal routines. CHFIV -- Cubic Hermite Function Integral Evaluator. (Real function called by PCHIA.) CHFMC -- Cubic Hermite Function Monotonicity Checker. (Integer function called by PCHMC.) PCHCE -- PCHIC End Derivative Setter. (Called by PCHIC.) PCHCI -- PCHIC Initial Derivative Setter. (Called by PCHIC.) PCHCS -- PCHIC Monotonicity Switch Derivative Setter. (Called by PCHIC.) PCHDF -- PCHIP Finite Difference Formula. (Real function called by PCHCE and PCHSP.) PCHST -- PCHIP Sign Testing Routine. (Real function called by various PCHIP routines.) PCHSW -- PCHCS Switch Excursion Adjuster. (Called by PCHCS.) The calling sequences for these routines are described in the prologues of the respective routines. To facilitate two-dimensional applications, the representation of a PCH function throughout the package includes INCFD, the in- crement between successive elements in the F- and D-arrays. For "normal" usage INCFD=1, and F and D are one-dimensional arrays. one would call PCHxx (where "xx" is "IM", "IC", or "SP") with N, X, F, D, 1 . Suppose, however, that one has data on a rectangular mesh, F2D(I,J) = value at (X(I), Y(J)), I=1(1)NX, J=1(1)NY. Assume the following dimensions: REAL X(NXMAX), Y(NYMAX) REAL F2D(NXMAX,NYMAX), FX(NXMAX,NYMAX), FY(NXMAX,NYMAX) where 2.LE.NX.LE.NXMAX AND 2.LE.NY.LE.NYMAX . To interpolate in X along the line Y = Y(J), call PCHxx with NX, X, F2D(1,J), FX(1,J), 1 . To interpolate along the line X = X(I), call PCHxx with NY, Y, F2D(I,1), FY(I,1), MXMAX . (This example assumes the usual columnwise storage of 2-D arrays in Fortran.) References [1] F.N.Fritsch and R.E.Carlson, "Monotone Piecewise Cubic Inter- polation," SIAM J. Numer. Anal. 17, 2 (April 1980), 238-246. [2] F.N.Fritsch and J.Butland, "A Method for Constructing Local Monotone Piecewise Cubic Interpolants," LLNL report UCRL-87559 (April 1982). [Submitted to Mathematics of Computation.] [3] F.N.Fritsch, "Piecewise Cubic Hermite Interpolation Package," LLNL report UCRL-87285 (July 1982). [Poster presented at the SIAM 30th Anniversary Meeting, 19-23 July 1982.] [4] Carl de Boor, A Practical Guide to Splines, Springer-Verlag (New York, 1978). [esp. Chapter IV, pp. 49-62.]
PCHFD
SUBROUTINE PCHFD(N,X,F,D,INCFD,SKIP,NE,XE,FE,DE,IERR) PCHFD: Piecewise Cubic Hermite Function and Derivative evaluator Evaluates the cubic Hermite function defined by N, X, F, D, to- gether with its first derivative, at the points XE(J), J=1(1)NE. If only function values are required, use PCHFE, instead. To provide compatibility with PCHIM and PCHIC, includes an increment between successive values of the F- and D-arrays.
Calling sequence: PARAMETER (INCFD = ...) INTEGER N, NE, IERR REAL X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE), DE(NE) LOGICAL SKIP CALL PCHFD (N, X, F, D, INCFD, SKIP, NE, XE, FE, DE, IERR) Parameters: N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real array of function values. F(1+(I-1)*INCFD) is the value corresponding to X(I). D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is the value corresponding to X(I). INCFD -- (input) increment between successive values in F and D. (Error return if INCFD.LT.1 .) SKIP -- (input/output) logical variable which should be set to .TRUE. if the user wishes to skip checks for validity of preceding parameters, or to .FALSE. otherwise. This will save time in case these checks have already been performed (say, in PCHIM or PCHIC). SKIP will be set to .TRUE. on normal return. NE -- (input) number of evaluation points. (Error return if NE.LT.1 .) XE -- (input) real array of points at which the functions are to be evaluated. NOTES: 1. The evaluation will be most efficient if the elements of XE are increasing relative to X; that is, XE(J) .GE. X(I) implies XE(K) .GE. X(I), all K.GE.J . 2. If any of the XE are outside the interval [X(1),X(N)], values are extrapolated from the nearest extreme cubic, and a warning error is returned. FE -- (output) real array of values of the cubic Hermite function defined by N, X, F, D at the points XE. DE -- (output) real array of values of the first derivative of the same function at the points XE. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). Warning error: IERR.GT.0 means that extrapolation was performed at IERR points. "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. IERR = -4 if NE.LT.1 . (Output arrays have not been changed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated. IERR = -5 if an error has occurred in the lower-level routine CHFDV. NB: this should never happen. Notify the author **IMMEDIATELY** if it does.PCHFE
SUBROUTINE PCHFE(N,X,F,D,INCFD,SKIP,NE,XE,FE,IERR) PCHFE: Piecewise Cubic Hermite Function Evaluator Evaluates the cubic Hermite function defined by N, X, F, D at the points XE(J), J=1(1)NE. To provide compatibility with PCHIM and PCHIC, includes an increment between successive values of the F- and D-arrays.
Calling sequence: PARAMETER (INCFD = ...) INTEGER N, NE, IERR REAL X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE) LOGICAL SKIP CALL PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR) Parameters: N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real array of function values. F(1+(I-1)*INCFD) is the value corresponding to X(I). D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is the value corresponding to X(I). INCFD -- (input) increment between successive values in F and D. (Error return if INCFD.LT.1 .) SKIP -- (input/output) logical variable which should be set to .TRUE. if the user wishes to skip checks for validity of preceding parameters, or to .FALSE. otherwise. This will save time in case these checks have already been performed (say, in PCHIM or PCHIC). SKIP will be set to .TRUE. on normal return. NE -- (input) number of evaluation points. (Error return if NE.LT.1 .) XE -- (input) real array of points at which the function is to be evaluated. NOTES: 1. The evaluation will be most efficient if the elements of XE are increasing relative to X; that is, XE(J) .GE. X(I) implies XE(K) .GE. X(I), all K.GE.J . 2. If any of the XE are outside the interval [X(1),X(N)], values are extrapolated from the nearest extreme cubic, and a warning error is returned. FE -- (output) real array of values of the cubic Hermite function defined by N, X, F, D at the points XE. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). Warning error: IERR.GT.0 means that extrapolation was performed at IERR points. "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. IERR = -4 if NE.LT.1 . (The FE-array has not been changed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated.PCHIA
REAL FUNCTION PCHIA(N,X,F,D,INCFD,SKIP,A,B,IERR) PCHIA: Piecewise Cubic Hermite Integrator, Arbitrary limits Evaluates the definite integral of the cubic Hermite function defined by N, X, F, D over the interval [A, B]. To provide compatibility with PCHIM and PCHIC, includes an increment between successive values of the F- and D-arrays.
Calling sequence: PARAMETER (INCFD = ...) INTEGER N, IERR REAL X(N), F(INCFD,N), D(INCFD,N), A, B LOGICAL SKIP VALUE = PCHIA (N, X, F, D, INCFD, SKIP, A, B, IERR) Parameters: VALUE -- (output) VALUE of the requested integral. N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real array of function values. F(1+(I-1)*INCFD) is the value corresponding to X(I). D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is the value corresponding to X(I). INCFD -- (input) increment between successive values in F and D. (Error return if INCFD.LT.1 .) SKIP -- (input/output) logical variable which should be set to .TRUE. if the user wishes to skip checks for validity of preceding parameters, or to .FALSE. otherwise. This will save time in case these checks have already been performed (say, in PCHIM or PCHIC). SKIP will be set to .TRUE. on return with IERR.GE.0 . A,B -- (input) the limits of integration. NOTE: There is no requirement that [A,B] be contained in [X(1),X(N)]. However, the resulting integral value will be highly suspect, if not. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). Warning errors: IERR = 1 if A is outside the interval [X(1),X(N)]. IERR = 2 if B is outside the interval [X(1),X(N)]. IERR = 3 if both of the above are true. (Note that this means that either [A,B] contains data interval or the intervals do not intersect at all.) "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. (Value has not been computed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated.PCHIC
SUBROUTINE PCHIC(IC,VC,SWITCH,N,X,F,D,INCFD,WK,NWK,IERR) PCHIC: Piecewise Cubic Hermite Interpolation Coefficients. Sets derivatives needed to determine a piecewise monotone piece- wise cubic interpolant to the data given in X and F satisfying the boundary conditions specified by IC and VC. The treatment of points where monotonicity switches direction is controlled by argument SWITCH. To facilitate two-dimensional applications, includes an increment between successive values of the F- and D-arrays. The resulting piecewise cubic Hermite function may be evaluated by PCHFE or PCHFD.
Calling sequence: PARAMETER (INCFD = ...) INTEGER IC(2), N, NWK, IERR REAL VC(2), SWITCH, X(N), F(INCFD,N), D(INCFD,N), WK(NWK) CALL PCHIC (IC, VC, SWITCH, N, X, F, D, INCFD, WK, NWK, IERR) Parameters: IC -- (input) integer array of length 2 specifying desired boundary conditions: IC(1) = IBEG, desired condition at beginning of data. IC(2) = IEND, desired condition at end of data. IBEG = 0 for the default boundary condition (the same as used by PCHIM). If IBEG.NE.0, then its sign indicates whether the boundary derivative is to be adjusted, if necessary, to be compatible with monotonicity: IBEG.GT.0 if no adjustment is to be performed. IBEG.LT.0 if the derivative is to be adjusted for monotonicity. Allowable values for the magnitude of IBEG are: IBEG = 1 if first derivative at X(1) is given in VC(1). IBEG = 2 if second derivative at X(1) is given in VC(1). IBEG = 3 to use the 3-point difference formula for D(1). (Reverts to the default b.c. if N.LT.3 .) IBEG = 4 to use the 4-point difference formula for D(1). (Reverts to the default b.c. if N.LT.4 .) IBEG = 5 to set D(1) so that the second derivative is con- tinuous at X(2). (Reverts to the default b.c. if N.LT.4.) This option is somewhat analogous to the "not a knot" boundary condition provided by PCHSP. NOTES (IBEG): 1. An error return is taken if ABS(IBEG).GT.5 . 2. Only in case IBEG.LE.0 is it guaranteed that the interpolant will be monotonic in the first interval. If the returned value of D(1) lies between zero and 3*SLOPE(1), the interpolant will be monotonic. This is **NOT** checked if IBEG.GT.0 . 3. If IBEG.LT.0 and D(1) had to be changed to achieve mono- tonicity, a warning error is returned. IEND may take on the same values as IBEG, but applied to derivative at X(N). In case IEND = 1 or 2, the value is given in VC(2). NOTES (IEND): 1. An error return is taken if ABS(IEND).GT.5 . 2. Only in case IEND.LE.0 is it guaranteed that the interpolant will be monotonic in the last interval. If the returned value of D(1+(N-1)*INCFD) lies between zero and 3*SLOPE(N-1), the interpolant will be monotonic. This is **NOT** checked if IEND.GT.0 . 3. If IEND.LT.0 and D(1+(N-1)*INCFD) had to be changed to achieve monotonicity, a warning error is returned. VC -- (input) real array of length 2 specifying desired boundary values, as indicated above. VC(1) need be set only if IC(1) = 1 or 2 . VC(2) need be set only if IC(2) = 1 or 2 . SWITCH -- (input) indicates desired treatment of points where direction of monotonicity switches: Set SWITCH to zero if interpolant is required to be mono- tonic in each interval, regardless of monotonicity of data. NOTES: 1. This will cause D to be set to zero at all switch points, thus forcing extrema there. 2. The result of using this option with the default boun- dary conditions will be identical to using PCHIM, but will generally cost more compute time. This option is provided only to facilitate comparison of different switch and/or boundary conditions. Set SWITCH nonzero to use a formula based on the 3-point difference formula in the vicinity of switch points. If SWITCH is positive, the interpolant on each interval containing an extremum is controlled to not deviate from the data by more than SWITCH*DFLOC, where DFLOC is the maximum of the change of F on this interval and its two immediate neighbors. If SWITCH is negative, no such control is to be imposed. N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real array of dependent variable values to be inter- polated. F(1+(I-1)*INCFD) is value corresponding to X(I). D -- (output) real array of derivative values at the data points. These values will determine a monotone cubic Hermite func- tion on each subinterval on which the data are monotonic, except possibly adjacent to switches in monotonicity. The value corresponding to X(I) is stored in D(1+(I-1)*INCFD), I=1(1)N. No other entries in D are changed. INCFD -- (input) increment between successive values in F and D. This argument is provided primarily for 2-D applications. (Error return if INCFD.LT.1 .) WK -- (scratch) real array of working storage. The user may wish to know that the returned values are: WK(I) = H(I) = X(I+1) - X(I) ; WK(N-1+I) = SLOPE(I) = (F(1,I+1) - F(1,I)) / H(I) for I = 1(1)N-1. NWK -- (input) length of work array. (Error return if NWK.LT.2*(N-1) .) IERR -- (output) error flag. Normal return: IERR = 0 (no errors). Warning errors: IERR = 1 if IBEG.LT.0 and D(1) had to be adjusted for monotonicity. IERR = 2 if IEND.LT.0 and D(1+(N-1)*INCFD) had to be adjusted for monotonicity. IERR = 3 if both of the above are true. "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. IERR = -4 if ABS(IBEG).GT.5 . IERR = -5 if ABS(IEND).GT.5 . IERR = -6 if both of the above are true. IERR = -7 if NWK.LT.2*(N-1) . (The D-array has not been changed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated.PCHID
REAL FUNCTION PCHID(N,X,F,D,INCFD,SKIP,IA,IB,IERR) PCHID: Piecewise Cubic Hermite Integrator, Data limits Evaluates the definite integral of the cubic Hermite function defined by N, X, F, D over the interval [X(IA), X(IB)]. To provide compatibility with PCHIM and PCHIC, includes an increment between successive values of the F- and D-arrays.
Calling sequence: PARAMETER (INCFD = ...) INTEGER N, IA, IB, IERR REAL X(N), F(INCFD,N), D(INCFD,N) LOGICAL SKIP VALUE = PCHID (N, X, F, D, INCFD, SKIP, IA, IB, IERR) Parameters: VALUE -- (output) VALUE of the requested integral. N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real array of function values. F(1+(I-1)*INCFD) is the value corresponding to X(I). D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is the value corresponding to X(I). INCFD -- (input) increment between successive values in F and D. (Error return if INCFD.LT.1 .) SKIP -- (input/output) logical variable which should be set to .TRUE. if the user wishes to skip checks for validity of preceding parameters, or to .FALSE. otherwise. This will save time in case these checks have already been performed (say, in PCHIM or PCHIC). SKIP will be set to .TRUE. on return with IERR = 0 or -4. IA,IB -- (input) indices in X-array for the limits of integration. both must be in the range [1,N]. (Error return if not.) No restrictions on their relative values. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. IERR = -4 if IA or IB is out of range. (Value has not been computed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated.PCHIM
SUBROUTINE PCHIM(N,X,F,D,INCFD,IERR) PCHIM: Piecewise Cubic Hermite Interpolation to Monotone data. Sets derivatives needed to determine a monotone piecewise cubic Hermite interpolant to the data given in X and F. Default boundary conditions are provided which are compatible with monotonicity. (See PCHIC if user control of boundary con- ditions is desired.) If the data are only piecewise monotonic, the interpolant will have an extremum at each point where monotonicity switches direc- tion. (See PCHIC if user control is desired in such cases.) To facilitate two-dimensional applications, includes an increment between successive values of the F- and D-arrays. The resulting piecewise cubic Hermite function may be evaluated by PCHFE or PCHFD.
Calling sequence: PARAMETER (INCFD = ...) INTEGER N, IERR REAL X(N), F(INCFD,N), D(INCFD,N) CALL PCHIM (N, X, F, D, INCFD, IERR) Parameters: N -- (input) number of data points. (Error return if N.LT.2 .) If N=2, simply does linear interpolation. X -- (input) real array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real array of dependent variable values to be inter- polated. F(1+(I-1)*INCFD) is value corresponding to X(I). PCHIM is designed for monotonic data, but it will work for any F-array. It will force extrema at points where mono- tonicity switches direction. If some other treatment of switch points is desired, PCHIC should be used instead. ----- D -- (output) real array of derivative values at the data points. If the data are monotonic, these values will determine a a monotone cubic Hermite function. The value corresponding to X(I) is stored in D(1+(I-1)*INCFD), I=1(1)N. No other entries in D are changed. INCFD -- (input) increment between successive values in F and D. This argument is provided primarily for 2-D applications. (Error return if INCFD.LT.1 .) IERR -- (output) error flag. Normal return: IERR = 0 (no errors). Warning error: IERR.GT.0 means that IERR switches in the direction of monotonicity were detected. "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. (The D-array has not been changed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated.PCHMC
SUBROUTINE PCHMC(N,X,F,D,INCFD,SKIP,ISMON,IERR) PCHMC: Piecewise Cubic Hermite Monotonicity Checker. Checks the cubic Hermite function defined by N, X, F, D for monotonicity. To provide compatibility with PCHIM and PCHIC, includes an increment between successive values of the F- and D-arrays.
Calling sequence: PARAMETER (INCFD = ...) INTEGER N, ISMON(N), IERR REAL X(N), F(INCFD,N), D(INCFD,N) LOGICAL SKIP CALL PCHMC (N, X, F, D, INCFD, SKIP, ISMON, IERR) Parameters: N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real array of function values. F(1+(I-1)*INCFD) is the value corresponding to X(I). D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is the value corresponding to X(I). INCFD -- (input) increment between successive values in F and D. (Error return if INCFD.LT.1 .) SKIP -- (input/output) logical variable which should be set to .TRUE. if the user wishes to skip checks for validity of preceding parameters, or to .FALSE. otherwise. This will save time in case these checks have already been performed. SKIP will be set to .TRUE. on normal return. ISMON -- (output) integer array indicating on which intervals the PCH function defined by N, X, F, D is monotonic. For data interval [X(I),X(I+1)], ISMON(I) = -1 if function is strictly decreasing; ISMON(I) = 0 if function is constant; ISMON(I) = 1 if function is strictly increasing; ISMON(I) = 2 if function is non-monotonic; ISMON(I) = 3 if unable to determine. (This means that the D-values are near the boundary of the monotonicity region. A small increase pro- duces non-monotonicity; decrease, strict monotonicity.) The above applies to I=1(1)N-1. ISMON(N) indicates whether the entire function is monotonic on [X(1),X(N)]. IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. (The ISMON-array has not been changed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated.PCHNGS
SUBROUTINE PCHNGS(II,XVAL,IPLACE,SX,IX,IRCX) THIS SUBROUTINE IS A MODIFICATION OF THE SUBROUTINE LCHNGS, SANDIA LABS. REPT. SAND78-0785. PCHNGS LIMITS THE TYPE OF STORAGE TO A SEQUENTIAL SCHEME. MODIFICATIONS BY K.L. HIEBERT AND R.J. HANSON REVISED 811130-1000 REVISED YYMMDD-HHMM SPARSE MATRIX ELEMENT ALTERATION SUBROUTINE. AUTHORS JOHN A. WISNIEWSKI AND R.J. HANSON SUBROUTINE PCHNGS() CHANGES ELEMENT II IN VECTOR +/- IRCX TO THE VALUE XVAL. II THE ABSOLUTE VALUE OF THIS INTEGER IS THE SUBSCRIPT FOR THE ELEMENT TO BE CHANGED. XVAL NEW VALUE OF THE MATRIX ELEMENT BEING CHANGED. IPLACE POINTER INFORMATION WHICH IS MAINTAINED BY THE PACKAGE. SX(*),IX(*) THE WORK ARRAYS WHICH ARE USED TO STORE THE SPARSE MATRIX. THESE ARRAYS ARE AUTOMATICALLY MAINTAINED BY THE PACKAGE FOR THE USER. IRCX POINTS TO THE VECTOR OF THE MATRIX BEING UPDATED. A NEGATIVE VALUE OF IRCX INDICATES THAT ROW -IRCX IS BEING UPDATED. A POSITIVE VALUE OF IRCX INDICATES THAT COLUMN IRCX IS BEING UPDATED. A ZERO VALUE OF IRCX IS AN ERROR. SINCE DATA ITEMS ARE KEPT SORTED IN THE SEQUENTIAL DATA STRUCTURE, CHANGING A MATRIX ELEMENT CAN REQUIRE THE MOVEMENT OF ALL THE DATA ITEMS IN THE MATRIX. FOR THIS REASON, IT IS SUGGESTED THAT DATA ITEMS BE ADDED A COL. AT A TIME, IN ASCENDING COL. SEQUENCE. FURTHERMORE, SINCE DELETING ITEMS FROM THE DATA STRUCTURE MAY ALSO REQUIRE MOVING LARGE AMOUNTS OF DATA, ZERO ELEMENTS ARE EXPLICITLY STORED IN THE MATRIX.
PCHSP
SUBROUTINE PCHSP(IC,VC,N,X,F,D,INCFD,WK,NWK,IERR) PCHSP: Piecewise Cubic Hermite Spline Computes the Hermite representation of the cubic spline inter- polant to the data given in X and F satisfying the boundary conditions specified by IC and VC. To facilitate two-dimensional applications, includes an increment between successive values of the F- and D-arrays. The resulting piecewise cubic Hermite function may be evaluated by PCHFE or PCHFD. NOTE: This is a modified version of C. de Boor'S cubic spline routine CUBSPL.
Calling sequence: PARAMETER (INCFD = ...) INTEGER IC(2), N, NWK, IERR REAL VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK) CALL PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR) Parameters: IC -- (input) integer array of length 2 specifying desired boundary conditions: IC(1) = IBEG, desired condition at beginning of data. IC(2) = IEND, desired condition at end of data. IBEG = 0 to set D(1) so that the third derivative is con- tinuous at X(2). This is the "not a knot" condition provided by de Boor'S cubic spline routine CUBSPL. < This is the default boundary condition. > IBEG = 1 if first derivative at X(1) is given in VC(1). IBEG = 2 if second derivative at X(1) is given in VC(1). IBEG = 3 to use the 3-point difference formula for D(1). (Reverts to the default b.c. if N.LT.3 .) IBEG = 4 to use the 4-point difference formula for D(1). (Reverts to the default b.c. if N.LT.4 .) NOTES: 1. An error return is taken if IBEG is out of range. 2. For the "natural" boundary condition, use IBEG=2 and VC(1)=0. IEND may take on the same values as IBEG, but applied to derivative at X(N). In case IEND = 1 or 2, the value is given in VC(2). NOTES: 1. An error return is taken if IEND is out of range. 2. For the "natural" boundary condition, use IEND=2 and VC(2)=0. VC -- (input) real array of length 2 specifying desired boundary values, as indicated above. VC(1) need be set only if IC(1) = 1 or 2 . VC(2) need be set only if IC(2) = 1 or 2 . N -- (input) number of data points. (Error return if N.LT.2 .) X -- (input) real array of independent variable values. The elements of X must be strictly increasing: X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) F -- (input) real array of dependent variable values to be inter- polated. F(1+(I-1)*INCFD) is value corresponding to X(I). D -- (output) real array of derivative values at the data points. These values will determine the cubic spline interpolant with the requested boundary conditions. The value corresponding to X(I) is stored in D(1+(I-1)*INCFD), I=1(1)N. No other entries in D are changed. INCFD -- (input) increment between successive values in F and D. This argument is provided primarily for 2-D applications. (Error return if INCFD.LT.1 .) WK -- (scratch) real array of working storage. NWK -- (input) length of work array. (Error return if NWK.LT.2*N .) IERR -- (output) error flag. Normal return: IERR = 0 (no errors). "Recoverable" errors: IERR = -1 if N.LT.2 . IERR = -2 if INCFD.LT.1 . IERR = -3 if the X-array is not strictly increasing. IERR = -4 if IBEG.LT.0 or IBEG.GT.4 . IERR = -5 if IEND.LT.0 of IEND.GT.4 . IERR = -6 if both of the above are true. IERR = -7 if NWK is too small. NOTE: The above errors are checked in the order listed, and following arguments have **NOT** been validated. (The D-array has not been changed in any of these cases.) IERR = -8 in case of trouble solving the linear system for the interior derivative values. (The D-array may have been changed in this case.) ( Do **NOT** use it! )PCHST
REAL FUNCTION PCHST(ARG1,ARG2) PCHST: PCHIP Sign-Testing Routine. Returns: -1. if ARG1 and ARG2 are of opposite sign. 0. if either argument is zero. +1. if ARG1 and ARG2 are of the same sign. The object is to do this without multiplying ARG1*ARG2, to avoid possible over/underflow problems. Fortran intrinsics used: SIGN.
Programmed by: Fred N. Fritsch, FTS 532-4275, (415) 422-4275, Mathematics and Statistics Division, Lawrence Livermore National Laboratory. Change record: 82-08-05 Converted to SLATEC library version.
Programming notes: To produce a double precision version, simply: a. Change PCHST to DPCHST wherever it occurs, b. Change all references to the Fortran intrinsics to their double presision equivalents, c. Change the real declarations to double precision, and d. Change the constants ZERO and ONE to double precision.PCHSW
SUBROUTINE PCHSW(DFMAX,IEXTRM,D1,D2,H,SLOPE,IERR) PCHSW: PCHCS Switch Excursion Limiter. Called by PCHCS to adjust D1 and D2 if necessary to insure that the extremum on this interval is not further than DFMAX from the extreme data value.
Calling sequence: INTEGER IEXTRM, IERR REAL DFMAX, D1, D2, H, SLOPE CALL PCHSW (DFMAX, IEXTRM, D1, D2, H, SLOPE, IERR) Parameters: DFMAX -- (input) maximum allowed difference between F(IEXTRM) and the cubic determined by derivative values D1,D2. (assumes DFMAX.GT.0.) IEXTRM -- (input) index of the extreme data value. (assumes IEXTRM = 1 or 2 . Any value .NE.1 is treated as 2.) D1,D2 -- (input) derivative values at the ends of the interval. (Assumes D1*D2 .LE. 0.) (output) may be modified if necessary to meet the restriction imposed by DFMAX. H -- (input) interval length. (Assumes H.GT.0.) SLOPE -- (input) data SLOPE on the interval. IERR -- (output) error flag. should be zero. If IERR=-1, assumption on D1 and D2 is not satisfied. If IERR=-2, quadratic equation locating extremum has negative descriminant (should never occur). ------- WARNING: This routine does no validity-checking of arguments. ------- Fortran intrinsics used: ABS, SIGN, SQRT.POLINT
SUBROUTINE POLINT(N,X,Y,C) Written by Robert E. Huddleston, Sandia Laboratories, Livermore Abstract Subroutine POLINT is designed to produce the polynomial which interpolates the data (X(I),Y(I)), I=1,...,N. POLINT sets up information in the array C which can be used by subroutine POLYVL to evaluate the polynomial and its derivatives and by subroutine POLCOF to produce the coefficients. Formal Parameters N - the number of data points (N .GE. 1) X - the array of abscissas (all of which must be distinct) Y - the array of ordinates C - an array of information used by subroutines ******* Dimensioning Information ******* Arrays X,Y, and C must be dimensioned at least N in the calling program.
POLYVL
SUBROUTINE POLYVL(NDER,XX,YFIT,YP,N,X,C,WORK,IERR) Written by Robert E. Huddleston, Sandia Laboratories, Livermore Abstract - Subroutine POLYVL calculates the value of the polynomial and its first NDER derivatives where the polynomial was produced by a previous call to HRMITE or POLINT. The variable N and the arrays X and C must not be altered between the call to HRMITE or POLINT and the call to POLYVL. ****** Dimensioning Information ******* YP must be dimensioned by at least NDER X must be dimensioned by at least N (see the abstract ) C must be dimensioned by at least N (see the abstract ) WORK must be dimensioned by at least 2*N if NDER is .GT. 0. *** Note *** If NDER=0, neither YP nor WORK need to be dimensioned variables. If NDER=1, YP does not need to be a dimensioned variable. ***** Input parameters NDER - the number of derivatives to be evaluated XX - the argument at which the polynomial and its derivatives are to be evaluated. N - ***** * N, X, and C must not be altered between the call X - * to HRMITE (or the call to POLINT) and the call * to POLYVL. C - ***** ***** Output Parameters YFIT - the value of the polynomial at XX YP - the derivatives of the polynomial at XX. The derivative of order J at XX is stored in YP(J) , J = 1,...,NDER. IERR - Output error flag with the following possible values. = 1 indicates normal execution ***** Storage Parameters WORK = this is an array to provide internal working storage for POLYVL. It must be dimensioned by at least 2*N if NDER is .GT. 0. If NDER=0, WORK does not need to be a dimensioned variable.
POS3D1
SUBROUTINE POS3D1(LP,L,MP,M,N,A,B,C,LDIMF,MDIMF,F,XRT,YRT,T,D,WX,
POSTG2
SUBROUTINE POSTG2(NPEROD,N,M,A,BB,C,IDIMQ,Q,B,B2,B3,W,W2,W3,D,
PPADD
SUBROUTINE PPADD(N,IERROR,A,C,CBP,BP,BH)
PPGQ8
SUBROUTINE PPGQ8(FUN,LDC,C,XI,LXI,KK,ID,A,B,INPPV,ERR,ANS,IERR)
PPGSF
FUNCTION PPGSF(X,IZ,C,A,BH)
PPPSF
FUNCTION PPPSF(X,IZ,C,A,BH)
PPVAL
FUNCTION PPVAL(LDC,C,XI,LXI,K,IDERIV,X,INPPV) Written by Carl de Boor and modified by D. E. Amos Reference SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. Abstract PPVAL is the PPVALU function of the reference. PPVAL calculates (at X) the value of the IDERIV-th derivative of the B-spline from the PP-representation (C,XI,LXI,K). The Taylor expansion about XI(J) for X in the interval XI(J) .LE. X .LT. XI(J+1) is evaluated, J=1,LXI. Right limiting values at X=XI(J) are obtained. PPVAL will extrapolate beyond XI(1) and XI(LXI+1). To obtain left limiting values (left derivatives) at XI(J), replace LXI by J-1 and set X=XI(J),J=2,LXI+1. PPVAL calls INTRV Description of Arguments Input LDC - leading dimension of C matrix, LDC .GE. K C - matrix of dimension at least (K,LXI) containing right derivatives at break points XI(*). XI - break point vector of length LXI+1 LXI - number of polynomial pieces K - order of B-spline, K .GE. 1 IDERIV - order of the derivative, 0 .LE. IDERIV .LE. K-1 IDERIV=0 gives the B-spline value X - argument, XI(1) .LE. X .LE. XI(LXI+1) INPPV - an initialization parameter which must be set to 1 the first time PPVAL is called. Output INPPV - INPPV contains information for efficient process- ing after the initial call and INPPV must not be changed by the user. Distinct splines require distinct INPPV parameters. PPVAL - value of the IDERIV-th derivative at X Error Conditions Improper input is a fatal error
PROC
SUBROUTINE PROC(ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,Y,M,A,B,C,D,W,U)
PROCP
SUBROUTINE PROCP(ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,Y,M,A,B,C,D,U,W)
PROD
SUBROUTINE PROD(ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,Y,M,A,B,C,D,W,U)
PRODP
SUBROUTINE PRODP(ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,Y,M,A,B,C,D,U,W)
PRVEC
FUNCTION PRVEC(M,U,V)
PRWPGE
SUBROUTINE PRWPGE(KEY,IPAGE,LPG,SX,IX) THIS SUBROUTINE IS A MODIFICATION OF THE SUBROUTINE LRWPGE, SANDIA LABS. REPT. SAND78-0785. PRWPGE LIMITS THE TYPE OF STORAGE TO A SEQUENTIAL SCHEME. MODIFICATIONS BY K.L. HIEBERT AND R.J. HANSON REVISED 811130-1000 REVISED YYMMDD-HHMM VIRTUAL MEMORY PAGE READ/WRITE SUBROUTINE. AUTHORS J. A. WISNIEWSKI AND R.J. HANSON DEPENDING ON THE VALUE OF KEY, SUBROUTINE PRWPGE() PERFORMS A PAGE READ OR WRITE OF PAGE IPAGE. THE PAGE HAS LENGTH LPG. KEY IS A FLAG INDICATING WHETHER A PAGE READ OR WRITE IS TO BE PERFORMED. IF KEY = 1 DATA IS READ. IF KEY = 2 DATA IS WRITTEN. IPAGE IS THE PAGE NUMBER OF THE MATRIX TO BE ACCESSED. LPG IS THE LENGTH OF THE PAGE OF THE MATRIX TO BE ACCESSED. SX(*),IX(*) IS THE MATRIX TO BE ACCESSED. EXTERNAL REFERENCES XERROR,PRWVIR
PRWVIR
SUBROUTINE PRWVIR(KEY,IPAGE,LPG,SX,IX) THIS SUBROUTINE IS A MODIFICATION OF THE SUBROUTINE LRWVIR. PRWVIR LIMITS THE TYPE OF STORAGE TO A SEQUENTIAL SPARSE MATRIX STORAGE SCHEME. THE PAGE STORAGE IS ON RANDOM ACCESS DISK. PRWVIR() IS PART OF THE SPARSE LP PACKAGE, SPLP( ). MODIFICATIONS BY K.L. HIEBERT AND R.J. HANSON AUTHORS J. A. WISNIEWSKI AND R.J. HANSON KEY IS A FLAG WHICH INDICATES WHETHER A READ OR WRITE OPERATION IS TO BE PERFORMED. A VALUE OF KEY=1 INDICATES A READ. A VALUE OF KEY=2 INDICATES A WRITE. IPAGE IS THE PAGE OF MATRIX MN WE ARE ACCESSING. LPG IS THE LENGTH OF THE PAGE. SX(*),IX(*) IS THE MATRIX DATA.
PSGF
FUNCTION PSGF(X,IZ,C,A,BH)