1*DECK DPCHBS 2 SUBROUTINE DPCHBS (N, X, F, D, INCFD, KNOTYP, NKNOTS, T, BCOEF, 3 + NDIM, KORD, IERR) 4C***BEGIN PROLOGUE DPCHBS 5C***PURPOSE Piecewise Cubic Hermite to B-Spline converter. 6C***LIBRARY SLATEC (PCHIP) 7C***CATEGORY E3 8C***TYPE DOUBLE PRECISION (PCHBS-S, DPCHBS-D) 9C***KEYWORDS B-SPLINES, CONVERSION, CUBIC HERMITE INTERPOLATION, 10C PIECEWISE CUBIC INTERPOLATION 11C***AUTHOR Fritsch, F. N., (LLNL) 12C Computing and Mathematics Research Division 13C Lawrence Livermore National Laboratory 14C P.O. Box 808 (L-316) 15C Livermore, CA 94550 16C FTS 532-4275, (510) 422-4275 17C***DESCRIPTION 18C 19C *Usage: 20C 21C INTEGER N, INCFD, KNOTYP, NKNOTS, NDIM, KORD, IERR 22C PARAMETER (INCFD = ...) 23C DOUBLE PRECISION X(nmax), F(INCFD,nmax), D(INCFD,nmax), 24C * T(2*nmax+4), BCOEF(2*nmax) 25C 26C CALL DPCHBS (N, X, F, D, INCFD, KNOTYP, NKNOTS, T, BCOEF, 27C * NDIM, KORD, IERR) 28C 29C *Arguments: 30C 31C N:IN is the number of data points, N.ge.2 . (not checked) 32C 33C X:IN is the real array of independent variable values. The 34C elements of X must be strictly increasing: 35C X(I-1) .LT. X(I), I = 2(1)N. (not checked) 36C nmax, the dimension of X, must be .ge.N. 37C 38C F:IN is the real array of dependent variable values. 39C F(1+(I-1)*INCFD) is the value corresponding to X(I). 40C nmax, the second dimension of F, must be .ge.N. 41C 42C D:IN is the real array of derivative values at the data points. 43C D(1+(I-1)*INCFD) is the value corresponding to X(I). 44C nmax, the second dimension of D, must be .ge.N. 45C 46C INCFD:IN is the increment between successive values in F and D. 47C This argument is provided primarily for 2-D applications. 48C It may have the value 1 for one-dimensional applications, 49C in which case F and D may be singly-subscripted arrays. 50C 51C KNOTYP:IN is a flag to control the knot sequence. 52C The knot sequence T is normally computed from X by putting 53C a double knot at each X and setting the end knot pairs 54C according to the value of KNOTYP: 55C KNOTYP = 0: Quadruple knots at X(1) and X(N). (default) 56C KNOTYP = 1: Replicate lengths of extreme subintervals: 57C T( 1 ) = T( 2 ) = X(1) - (X(2)-X(1)) ; 58C T(M+4) = T(M+3) = X(N) + (X(N)-X(N-1)). 59C KNOTYP = 2: Periodic placement of boundary knots: 60C T( 1 ) = T( 2 ) = X(1) - (X(N)-X(N-1)); 61C T(M+4) = T(M+3) = X(N) + (X(2)-X(1)) . 62C Here M=NDIM=2*N. 63C If the input value of KNOTYP is negative, however, it is 64C assumed that NKNOTS and T were set in a previous call. 65C This option is provided for improved efficiency when used 66C in a parametric setting. 67C 68C NKNOTS:INOUT is the number of knots. 69C If KNOTYP.GE.0, then NKNOTS will be set to NDIM+4. 70C If KNOTYP.LT.0, then NKNOTS is an input variable, and an 71C error return will be taken if it is not equal to NDIM+4. 72C 73C T:INOUT is the array of 2*N+4 knots for the B-representation. 74C If KNOTYP.GE.0, T will be returned by DPCHBS with the 75C interior double knots equal to the X-values and the 76C boundary knots set as indicated above. 77C If KNOTYP.LT.0, it is assumed that T was set by a 78C previous call to DPCHBS. (This routine does **not** 79C verify that T forms a legitimate knot sequence.) 80C 81C BCOEF:OUT is the array of 2*N B-spline coefficients. 82C 83C NDIM:OUT is the dimension of the B-spline space. (Set to 2*N.) 84C 85C KORD:OUT is the order of the B-spline. (Set to 4.) 86C 87C IERR:OUT is an error flag. 88C Normal return: 89C IERR = 0 (no errors). 90C "Recoverable" errors: 91C IERR = -4 if KNOTYP.GT.2 . 92C IERR = -5 if KNOTYP.LT.0 and NKNOTS.NE.(2*N+4). 93C 94C *Description: 95C DPCHBS computes the B-spline representation of the PCH function 96C determined by N,X,F,D. To be compatible with the rest of PCHIP, 97C DPCHBS includes INCFD, the increment between successive values of 98C the F- and D-arrays. 99C 100C The output is the B-representation for the function: NKNOTS, T, 101C BCOEF, NDIM, KORD. 102C 103C *Caution: 104C Since it is assumed that the input PCH function has been 105C computed by one of the other routines in the package PCHIP, 106C input arguments N, X, INCFD are **not** checked for validity. 107C 108C *Restrictions/assumptions: 109C 1. N.GE.2 . (not checked) 110C 2. X(i).LT.X(i+1), i=1,...,N . (not checked) 111C 3. INCFD.GT.0 . (not checked) 112C 4. KNOTYP.LE.2 . (error return if not) 113C *5. NKNOTS = NDIM+4 = 2*N+4 . (error return if not) 114C *6. T(2*k+1) = T(2*k) = X(k), k=1,...,N . (not checked) 115C 116C * Indicates this applies only if KNOTYP.LT.0 . 117C 118C *Portability: 119C Argument INCFD is used only to cause the compiler to generate 120C efficient code for the subscript expressions (1+(I-1)*INCFD) . 121C The normal usage, in which DPCHBS is called with one-dimensional 122C arrays F and D, is probably non-Fortran 77, in the strict sense, 123C but it works on all systems on which DPCHBS has been tested. 124C 125C *See Also: 126C PCHIC, PCHIM, or PCHSP can be used to determine an interpolating 127C PCH function from a set of data. 128C The B-spline routine DBVALU can be used to evaluate the 129C B-representation that is output by DPCHBS. 130C (See BSPDOC for more information.) 131C 132C***REFERENCES F. N. Fritsch, "Representations for parametric cubic 133C splines," Computer Aided Geometric Design 6 (1989), 134C pp.79-82. 135C***ROUTINES CALLED DPCHKT, XERMSG 136C***REVISION HISTORY (YYMMDD) 137C 870701 DATE WRITTEN 138C 900405 Converted Fortran to upper case. 139C 900405 Removed requirement that X be dimensioned N+1. 140C 900406 Modified to make PCHKT a subsidiary routine to simplify 141C usage. In the process, added argument INCFD to be com- 142C patible with the rest of PCHIP. 143C 900410 Converted prologue to SLATEC 4.0 format. 144C 900410 Added calls to XERMSG and changed constant 3. to 3 to 145C reduce single/double differences. 146C 900411 Added reference. 147C 900430 Produced double precision version. 148C 900501 Corrected declarations. 149C 930317 Minor cosmetic changes. (FNF) 150C 930514 Corrected problems with dimensioning of arguments and 151C clarified DESCRIPTION. (FNF) 152C 930604 Removed NKNOTS from DPCHKT call list. (FNF) 153C***END PROLOGUE DPCHBS 154C 155C*Internal Notes: 156C 157C**End 158C 159C Declare arguments. 160C 161 INTEGER N, INCFD, KNOTYP, NKNOTS, NDIM, KORD, IERR 162 DOUBLE PRECISION X(*), F(INCFD,*), D(INCFD,*), T(*), BCOEF(*) 163C 164C Declare local variables. 165C 166 INTEGER K, KK 167 DOUBLE PRECISION DOV3, HNEW, HOLD 168 CHARACTER*8 LIBNAM, SUBNAM 169C***FIRST EXECUTABLE STATEMENT DPCHBS 170C 171C Initialize. 172C 173 NDIM = 2*N 174 KORD = 4 175 IERR = 0 176 LIBNAM = 'SLATEC' 177 SUBNAM = 'DPCHBS' 178C 179C Check argument validity. Set up knot sequence if OK. 180C 181 IF ( KNOTYP.GT.2 ) THEN 182 IERR = -1 183 CALL XERMSG (LIBNAM, SUBNAM, 'KNOTYP GREATER THAN 2', IERR, 1) 184 RETURN 185 ENDIF 186 IF ( KNOTYP.LT.0 ) THEN 187 IF ( NKNOTS.NE.NDIM+4 ) THEN 188 IERR = -2 189 CALL XERMSG (LIBNAM, SUBNAM, 190 * 'KNOTYP.LT.0 AND NKNOTS.NE.(2*N+4)', IERR, 1) 191 RETURN 192 ENDIF 193 ELSE 194C Set up knot sequence. 195 NKNOTS = NDIM + 4 196 CALL DPCHKT (N, X, KNOTYP, T) 197 ENDIF 198C 199C Compute B-spline coefficients. 200C 201 HNEW = T(3) - T(1) 202 DO 40 K = 1, N 203 KK = 2*K 204 HOLD = HNEW 205C The following requires mixed mode arithmetic. 206 DOV3 = D(1,K)/3 207 BCOEF(KK-1) = F(1,K) - HOLD*DOV3 208C The following assumes T(2*K+1) = X(K). 209 HNEW = T(KK+3) - T(KK+1) 210 BCOEF(KK) = F(1,K) + HNEW*DOV3 211 40 CONTINUE 212C 213C Terminate. 214C 215 RETURN 216C------------- LAST LINE OF DPCHBS FOLLOWS ----------------------------- 217 END 218