*DECK PCHBS
SUBROUTINE PCHBS (N, X, F, D, INCFD, KNOTYP, NKNOTS, T, BCOEF,
+ NDIM, KORD, IERR)
C***BEGIN PROLOGUE PCHBS
C***PURPOSE Piecewise Cubic Hermite to B-Spline converter.
C***LIBRARY SLATEC (PCHIP)
C***CATEGORY E3
C***TYPE SINGLE PRECISION (PCHBS-S, DPCHBS-D)
C***KEYWORDS B-SPLINES, CONVERSION, CUBIC HERMITE INTERPOLATION,
C PIECEWISE CUBIC INTERPOLATION
C***AUTHOR Fritsch, F. N., (LLNL)
C Computing and Mathematics Research Division
C Lawrence Livermore National Laboratory
C P.O. Box 808 (L-316)
C Livermore, CA 94550
C FTS 532-4275, (510) 422-4275
C***DESCRIPTION
C
C *Usage:
C
C INTEGER N, INCFD, KNOTYP, NKNOTS, NDIM, KORD, IERR
C PARAMETER (INCFD = ...)
C REAL X(nmax), F(INCFD,nmax), D(INCFD,nmax), T(2*nmax+4),
C * BCOEF(2*nmax)
C
C CALL PCHBS (N, X, F, D, INCFD, KNOTYP, NKNOTS, T, BCOEF,
C * NDIM, KORD, IERR)
C
C *Arguments:
C
C N:IN is the number of data points, N.ge.2 . (not checked)
C
C X:IN is the real array of independent variable values. The
C elements of X must be strictly increasing:
C X(I-1) .LT. X(I), I = 2(1)N. (not checked)
C nmax, the dimension of X, must be .ge.N.
C
C F:IN is the real array of dependent variable values.
C F(1+(I-1)*INCFD) is the value corresponding to X(I).
C nmax, the second dimension of F, must be .ge.N.
C
C D:IN is the real array of derivative values at the data points.
C D(1+(I-1)*INCFD) is the value corresponding to X(I).
C nmax, the second dimension of D, must be .ge.N.
C
C INCFD:IN is the increment between successive values in F and D.
C This argument is provided primarily for 2-D applications.
C It may have the value 1 for one-dimensional applications,
C in which case F and D may be singly-subscripted arrays.
C
C KNOTYP:IN is a flag to control the knot sequence.
C The knot sequence T is normally computed from X by putting
C a double knot at each X and setting the end knot pairs
C according to the value of KNOTYP:
C KNOTYP = 0: Quadruple knots at X(1) and X(N). (default)
C KNOTYP = 1: Replicate lengths of extreme subintervals:
C T( 1 ) = T( 2 ) = X(1) - (X(2)-X(1)) ;
C T(M+4) = T(M+3) = X(N) + (X(N)-X(N-1)).
C KNOTYP = 2: Periodic placement of boundary knots:
C T( 1 ) = T( 2 ) = X(1) - (X(N)-X(N-1));
C T(M+4) = T(M+3) = X(N) + (X(2)-X(1)) .
C Here M=NDIM=2*N.
C If the input value of KNOTYP is negative, however, it is
C assumed that NKNOTS and T were set in a previous call.
C This option is provided for improved efficiency when used
C in a parametric setting.
C
C NKNOTS:INOUT is the number of knots.
C If KNOTYP.GE.0, then NKNOTS will be set to NDIM+4.
C If KNOTYP.LT.0, then NKNOTS is an input variable, and an
C error return will be taken if it is not equal to NDIM+4.
C
C T:INOUT is the array of 2*N+4 knots for the B-representation.
C If KNOTYP.GE.0, T will be returned by PCHBS with the
C interior double knots equal to the X-values and the
C boundary knots set as indicated above.
C If KNOTYP.LT.0, it is assumed that T was set by a
C previous call to PCHBS. (This routine does **not**
C verify that T forms a legitimate knot sequence.)
C
C BCOEF:OUT is the array of 2*N B-spline coefficients.
C
C NDIM:OUT is the dimension of the B-spline space. (Set to 2*N.)
C
C KORD:OUT is the order of the B-spline. (Set to 4.)
C
C IERR:OUT is an error flag.
C Normal return:
C IERR = 0 (no errors).
C "Recoverable" errors:
C IERR = -4 if KNOTYP.GT.2 .
C IERR = -5 if KNOTYP.LT.0 and NKNOTS.NE.(2*N+4).
C
C *Description:
C PCHBS computes the B-spline representation of the PCH function
C determined by N,X,F,D. To be compatible with the rest of PCHIP,
C PCHBS includes INCFD, the increment between successive values of
C the F- and D-arrays.
C
C The output is the B-representation for the function: NKNOTS, T,
C BCOEF, NDIM, KORD.
C
C *Caution:
C Since it is assumed that the input PCH function has been
C computed by one of the other routines in the package PCHIP,
C input arguments N, X, INCFD are **not** checked for validity.
C
C *Restrictions/assumptions:
C 1. N.GE.2 . (not checked)
C 2. X(i).LT.X(i+1), i=1,...,N . (not checked)
C 3. INCFD.GT.0 . (not checked)
C 4. KNOTYP.LE.2 . (error return if not)
C *5. NKNOTS = NDIM+4 = 2*N+4 . (error return if not)
C *6. T(2*k+1) = T(2*k) = X(k), k=1,...,N . (not checked)
C
C * Indicates this applies only if KNOTYP.LT.0 .
C
C *Portability:
C Argument INCFD is used only to cause the compiler to generate
C efficient code for the subscript expressions (1+(I-1)*INCFD) .
C The normal usage, in which PCHBS is called with one-dimensional
C arrays F and D, is probably non-Fortran 77, in the strict sense,
C but it works on all systems on which PCHBS has been tested.
C
C *See Also:
C PCHIC, PCHIM, or PCHSP can be used to determine an interpolating
C PCH function from a set of data.
C The B-spline routine BVALU can be used to evaluate the
C B-representation that is output by PCHBS.
C (See BSPDOC for more information.)
C
C***REFERENCES F. N. Fritsch, "Representations for parametric cubic
C splines," Computer Aided Geometric Design 6 (1989),
C pp.79-82.
C***ROUTINES CALLED PCHKT, XERMSG
C***REVISION HISTORY (YYMMDD)
C 870701 DATE WRITTEN
C 900405 Converted Fortran to upper case.
C 900405 Removed requirement that X be dimensioned N+1.
C 900406 Modified to make PCHKT a subsidiary routine to simplify
C usage. In the process, added argument INCFD to be com-
C patible with the rest of PCHIP.
C 900410 Converted prologue to SLATEC 4.0 format.
C 900410 Added calls to XERMSG and changed constant 3. to 3 to
C reduce single/double differences.
C 900411 Added reference.
C 900501 Corrected declarations.
C 930317 Minor cosmetic changes. (FNF)
C 930514 Corrected problems with dimensioning of arguments and
C clarified DESCRIPTION. (FNF)
C 930604 Removed NKNOTS from PCHKT call list. (FNF)
C***END PROLOGUE PCHBS
C
C*Internal Notes:
C
C**End
C
C Declare arguments.
C
INTEGER N, INCFD, KNOTYP, NKNOTS, NDIM, KORD, IERR
REAL X(*), F(INCFD,*), D(INCFD,*), T(*), BCOEF(*)
C
C Declare local variables.
C
INTEGER K, KK
REAL DOV3, HNEW, HOLD
CHARACTER*8 LIBNAM, SUBNAM
C***FIRST EXECUTABLE STATEMENT PCHBS
C
C Initialize.
C
NDIM = 2*N
KORD = 4
IERR = 0
LIBNAM = 'SLATEC'
SUBNAM = 'PCHBS'
C
C Check argument validity. Set up knot sequence if OK.
C
IF ( KNOTYP.GT.2 ) THEN
IERR = -1
CALL XERMSG (LIBNAM, SUBNAM, 'KNOTYP GREATER THAN 2', IERR, 1)
RETURN
ENDIF
IF ( KNOTYP.LT.0 ) THEN
IF ( NKNOTS.NE.NDIM+4 ) THEN
IERR = -2
CALL XERMSG (LIBNAM, SUBNAM,
* 'KNOTYP.LT.0 AND NKNOTS.NE.(2*N+4)', IERR, 1)
RETURN
ENDIF
ELSE
C Set up knot sequence.
NKNOTS = NDIM + 4
CALL PCHKT (N, X, KNOTYP, T)
ENDIF
C
C Compute B-spline coefficients.
C
HNEW = T(3) - T(1)
DO 40 K = 1, N
KK = 2*K
HOLD = HNEW
C The following requires mixed mode arithmetic.
DOV3 = D(1,K)/3
BCOEF(KK-1) = F(1,K) - HOLD*DOV3
C The following assumes T(2*K+1) = X(K).
HNEW = T(KK+3) - T(KK+1)
BCOEF(KK) = F(1,K) + HNEW*DOV3
40 CONTINUE
C
C Terminate.
C
RETURN
C------------- LAST LINE OF PCHBS FOLLOWS ------------------------------
END