!DECK PCHSP SUBROUTINE PCHSP(IC, VC, N, X, F, D, INCFD, WK, NWK, IERR) !***BEGIN PROLOGUE PCHSP !***PURPOSE Set derivatives needed to determine the Hermite represen- ! tation of the cubic spline interpolant to given data, with ! specified boundary conditions. !***LIBRARY SLATEC (PCHIP) !***CATEGORY E1A !***TYPE SINGLE PRECISION (PCHSP-S, DPCHSP-D) !***KEYWORDS CUBIC HERMITE INTERPOLATION, PCHIP, ! PIECEWISE CUBIC INTERPOLATION, SPLINE INTERPOLATION !***AUTHOR Fritsch, F. N., (LLNL) ! Lawrence Livermore National Laboratory ! P.O. Box 808 (L-316) ! Livermore, CA 94550 ! FTS 532-4275, (510) 422-4275 !***DESCRIPTION ! ! 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! ) ! !***REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- ! Verlag, New York, 1978, pp. 53-59. !***ROUTINES CALLED PCHDF, XERMSG !***REVISION HISTORY (YYMMDD) ! 820503 DATE WRITTEN ! 820804 Converted to SLATEC library version. ! 870707 Minor cosmetic changes to prologue. ! 890411 Added SAVE statements (Vers. 3.2). ! 890703 Corrected category record. (WRB) ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) ! 920429 Revised format and order of references. (WRB,FNF) !***END PROLOGUE PCHSP ! Programming notes: ! ! To produce a double precision version, simply: ! a. Change PCHSP to DPCHSP wherever it occurs, ! b. Change the real declarations to double precision, and ! c. Change the constants ZERO, HALF, ... to double precision. ! ! DECLARE ARGUMENTS. ! INTEGER :: IC(2), N, INCFD, NWK, IERR REAL :: VC(2), X(*), F(INCFD,*), D(INCFD,*), WK(2,*) ! ! DECLARE LOCAL VARIABLES. ! INTEGER :: IBEG, IEND, INDEX, J, NM1 REAL :: G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), ZERO SAVE ZERO, HALF, ONE, TWO, THREE REAL :: PCHDF ! DATA ZERO /0./, HALF /0.5/, ONE /1./, TWO /2./, THREE /3./ ! ! VALIDITY-CHECK ARGUMENTS. ! !***FIRST EXECUTABLE STATEMENT PCHSP IF ( N<2 ) GO TO 5001 IF ( INCFD<1 ) GO TO 5002 DO J = 2, N IF ( X(J)<=X(J-1) ) GO TO 5003 END DO ! IBEG = IC(1) IEND = IC(2) IERR = 0 IF ( (IBEG<0).OR.(IBEG>4) ) IERR = IERR - 1 IF ( (IEND<0).OR.(IEND>4) ) IERR = IERR - 2 IF ( IERR<0 ) GO TO 5004 ! ! FUNCTION DEFINITION IS OK -- GO ON. ! IF ( NWK < 2*N ) GO TO 5007 ! ! COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO, ! COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.). DO J=2,N WK(1,J) = X(J) - X(J-1) WK(2,J) = (F(1,J) - F(1,J-1))/WK(1,J) END DO ! ! SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. ! IF ( IBEG>N ) IBEG = 0 IF ( IEND>N ) IEND = 0 ! ! SET UP FOR BOUNDARY CONDITIONS. ! IF ( (IBEG==1).OR.(IBEG==2) ) THEN D(1,1) = VC(1) ELSE IF (IBEG > 2) THEN ! PICK UP FIRST IBEG POINTS, IN REVERSE ORDER. DO J = 1, IBEG INDEX = IBEG-J+1 ! INDEX RUNS FROM IBEG DOWN TO 1. XTEMP(J) = X(INDEX) IF (J < IBEG) STEMP(J) = WK(2,INDEX) END DO ! -------------------------------- D(1,1) = PCHDF (IBEG, XTEMP, STEMP, IERR) ! -------------------------------- IF (IERR /= 0) GO TO 5009 IBEG = 1 ENDIF ! IF ( (IEND==1).OR.(IEND==2) ) THEN D(1,N) = VC(2) ELSE IF (IEND > 2) THEN ! PICK UP LAST IEND POINTS. DO J = 1, IEND INDEX = N-IEND+J ! INDEX RUNS FROM N+1-IEND UP TO N. XTEMP(J) = X(INDEX) IF (J < IEND) STEMP(J) = WK(2,INDEX+1) END DO ! -------------------------------- D(1,N) = PCHDF (IEND, XTEMP, STEMP, IERR) ! -------------------------------- IF (IERR /= 0) GO TO 5009 IEND = 1 ENDIF ! ! --------------------( BEGIN CODING FROM CUBSPL )-------------------- ! ! **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF ! F AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM- ! INATION, WITH S(J) ENDING UP IN D(1,J), ALL J. ! WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE. ! ! CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM ! WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1) ! IF (IBEG == 0) THEN IF (N == 2) THEN ! NO CONDITION AT LEFT END AND N = 2. WK(2,1) = ONE WK(1,1) = ONE D(1,1) = TWO*WK(2,2) ELSE ! NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2. WK(2,1) = WK(1,3) WK(1,1) = WK(1,2) + WK(1,3) D(1,1) =((WK(1,2) + TWO*WK(1,1))*WK(2,2)*WK(1,3) & + WK(1,2)**2*WK(2,3)) / WK(1,1) ENDIF ELSE IF (IBEG == 1) THEN ! SLOPE PRESCRIBED AT LEFT END. WK(2,1) = ONE WK(1,1) = ZERO ELSE ! SECOND DERIVATIVE PRESCRIBED AT LEFT END. WK(2,1) = TWO WK(1,1) = ONE D(1,1) = THREE*WK(2,2) - HALF*WK(1,2)*D(1,1) ENDIF ! ! IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND ! CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH ! EQUATION READS WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J). ! NM1 = N-1 IF (NM1 > 1) THEN DO J=2,NM1 IF (WK(2,J-1) == ZERO) GO TO 5008 G = -WK(1,J+1)/WK(2,J-1) D(1,J) = G*D(1,J-1) & + THREE*(WK(1,J)*WK(2,J+1) + WK(1,J+1)*WK(2,J)) WK(2,J) = G*WK(1,J-1) + TWO*(WK(1,J) + WK(1,J+1)) END DO ENDIF ! ! CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM ! (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N) ! ! IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK- ! SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT ! AT THIS POINT. IF (IEND == 1) GO TO 30 ! IF (IEND == 0) THEN IF (N==2 .AND. IBEG==0) THEN ! NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2. D(1,2) = WK(2,2) GO TO 30 ELSE IF ((N==2) .OR. (N==3 .AND. IBEG==0)) THEN ! EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT* ! NOT-A-KNOT AT LEFT END POINT). D(1,N) = TWO*WK(2,N) WK(2,N) = ONE IF (WK(2,N-1) == ZERO) GO TO 5008 G = -ONE/WK(2,N-1) ELSE ! NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR ALSO NOT-A- ! KNOT AT LEFT END POINT. G = WK(1,N-1) + WK(1,N) ! DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES). D(1,N) = ((WK(1,N)+TWO*G)*WK(2,N)*WK(1,N-1) & + WK(1,N)**2*(F(1,N-1)-F(1,N-2))/WK(1,N-1))/G IF (WK(2,N-1) == ZERO) GO TO 5008 G = -G/WK(2,N-1) WK(2,N) = WK(1,N-1) ENDIF ELSE ! SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT. D(1,N) = THREE*WK(2,N) + HALF*WK(1,N)*D(1,N) WK(2,N) = TWO IF (WK(2,N-1) == ZERO) GO TO 5008 G = -ONE/WK(2,N-1) ENDIF ! ! COMPLETE FORWARD PASS OF GAUSS ELIMINATION. ! WK(2,N) = G*WK(1,N-1) + WK(2,N) IF (WK(2,N) == ZERO) GO TO 5008 D(1,N) = (G*D(1,N-1) + D(1,N))/WK(2,N) ! ! CARRY OUT BACK SUBSTITUTION ! 30 CONTINUE DO J=NM1,1,-1 IF (WK(2,J) == ZERO) GO TO 5008 D(1,J) = (D(1,J) - WK(1,J)*D(1,J+1))/WK(2,J) END DO ! --------------------( END CODING FROM CUBSPL )-------------------- ! ! NORMAL RETURN. ! RETURN ! ! ERROR RETURNS. ! 5001 CONTINUE ! N.LT.2 RETURN. IERR = -1 CALL XERMSG ('SLATEC', 'PCHSP', & 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1) RETURN ! 5002 CONTINUE ! INCFD.LT.1 RETURN. IERR = -2 CALL XERMSG ('SLATEC', 'PCHSP', 'INCREMENT LESS THAN ONE', IERR, & 1) RETURN ! 5003 CONTINUE ! X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERMSG ('SLATEC', 'PCHSP', 'X-ARRAY NOT STRICTLY INCREASING' & , IERR, 1) RETURN ! 5004 CONTINUE ! IC OUT OF RANGE RETURN. IERR = IERR - 3 CALL XERMSG ('SLATEC', 'PCHSP', 'IC OUT OF RANGE', IERR, 1) RETURN ! 5007 CONTINUE ! NWK TOO SMALL RETURN. IERR = -7 CALL XERMSG ('SLATEC', 'PCHSP', 'WORK ARRAY TOO SMALL', IERR, 1) RETURN ! 5008 CONTINUE ! SINGULAR SYSTEM. ! *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES *** ! *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). *** IERR = -8 CALL XERMSG ('SLATEC', 'PCHSP', 'SINGULAR LINEAR SYSTEM', IERR, & 1) RETURN ! 5009 CONTINUE ! ERROR RETURN FROM PCHDF. ! *** THIS CASE SHOULD NEVER OCCUR *** IERR = -9 CALL XERMSG ('SLATEC', 'PCHSP', 'ERROR RETURN FROM PCHDF', IERR, & 1) RETURN !------------- LAST LINE OF PCHSP FOLLOWS ------------------------------ END SUBROUTINE PCHSP