SUBROUTINE DPCHEZ(N,X,F,D,SPLINE,WK,LWK,IERR) C***BEGIN PROLOGUE DPCHEZ C***DATE WRITTEN 870821 (YYMMDD) C***REVISION DATE 870908 (YYMMDD) C***CATEGORY NO. E1B C***KEYWORDS CUBIC HERMITE MONOTONE INTERPOLATION, SPLINE C INTERPOLATION, EASY TO USE PIECEWISE CUBIC INTERPOLATION C***AUTHOR KAHANER, D.K., (NBS) C SCIENTIFIC COMPUTING DIVISION C NATIONAL BUREAU OF STANDARDS C GAITHERSBURG, MARYLAND 20899 C (301) 975-3808 C***PURPOSE Easy to use spline or cubic Hermite interpolation. C***DESCRIPTION C C DPCHEZ: Piecewise Cubic Interpolation, Easy to Use. C C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C C Sets derivatives for spline (two continuous derivatives) or C Hermite cubic (one continuous derivative) interpolation. C Spline interpolation is smoother, but may not "look" right if the C data contains both "steep" and "flat" sections. Hermite cubics C can produce a "visually pleasing" and monotone interpolant to C monotone data. This is an easy to use driver for the routines C by F. N. Fritsch in reference (4) below. Various boundary C conditions are set to default values by DPCHEZ. Many other choices C are available in the subroutines PCHIC, DPCHIM and DPCHSP. C C Use PCHEV to evaluate the resulting function and its derivative. C C ---------------------------------------------------------------------- C C Calling sequence: CALL DPCHEZ (N, X, F, D, SPLINE, WK, LWK, IERR) C C INTEGER N, IERR, LWK C DOUBLE PRECISION X(N), F(N), D(N), WK(*) C LOGICAL SPLINE C C Parameters: C C N -- (input) number of data points. (Error return if N.LT.2 .) C If N=2, simply does linear interpolation. C C X -- (input) 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. C (Error return if not.) C C F -- (input) real array of dependent variable values to be inter- C polated. F(I) is value corresponding to X(I). C C D -- (output) real array of derivative values at the data points. C C SPLINE -- (input) logical variable to specify if the interpolant C is to be a spline with two continuous derivaties C (set SPLINE=.TRUE.) or a Hermite cubic interpolant with one C continuous derivative (set SPLINE=.FALSE.). C Note: If SPLINE=.TRUE. the interpolating spline satisfies the C default "not-a-knot" boundary condition, with a continuous C third derivative at X(2) and X(N-1). See reference (3). C If SPLINE=.FALSE. the interpolating Hermite cubic will be C monotone if the input data is monotone. Boundary conditions C computed from the derivative of a local quadratic unless thi C alters monotonicity. C C WK -- (scratch) real work array, which must be declared by the cal C program to be at least 2*N if SPLINE is .TRUE. and not used C otherwise. C C LWK -- (input) length of work array WK. (Error return if C LWK.LT.2*N and SPLINE is .TRUE., not checked otherwise.) C C IERR -- (output) error flag. C Normal return: C IERR = 0 (no errors). C Warning error: C IERR.GT.0 (can only occur when SPLINE=.FALSE.) means tha C IERR switches in the direction of monotonicity were de C When SPLINE=.FALSE., DPCHEZ guarantees that if the inp C data is monotone, the interpolant will be too. This wa C is to alert you to the fact that the input data was no C monotone. C "Recoverable" errors: C IERR = -1 if N.LT.2 . C IERR = -3 if the X-array is not strictly increasing. C IERR = -7 if LWK is less than 2*N and SPLINE is .TRUE. C (The D-array has not been changed in any of these cases.) C NOTE: The above errors are checked in the order listed, C and following arguments have **NOT** been validated. C C ---------------------------------------------------------------------- C***REFERENCES 1. F.N.FRITSCH AND R.E.CARLSON, 'MONOTONE PIECEWISE C CUBIC INTERPOLATION,' SIAM J.NUMER.ANAL. 17, 2 (APRIL C 1980), 238-246. C 2. F.N.FRITSCH AND J.BUTLAND, 'A METHOD FOR CONSTRUCTING C LOCAL MONOTONE PIECEWISE CUBIC INTERPOLANTS,' LLNL C PREPRINT UCRL-87559 (APRIL 1982). C 3. CARL DE BOOR, A PRACTICAL GUIDE TO SPLINES, SPRINGER- C VERLAG (NEW YORK, 1978). (ESP. CHAPTER IV, PP.49-62.) C 4. F.N.FRITSCH, 'PIECEWISE CUBIC HERMITE INTERPOLATION C PACKAGE, FINAL SPECIFICATIONS', LAWRENCE LIVERMORE C NATIONAL LABORATORY, COMPUTER DOCUMENTATION UCID-30194 C AUGUST 1982. C***ROUTINES CALLED DPCHIM,DPCHSP C***END PROLOGUE DPCHEZ INTEGER N, LWK, IERR DOUBLE PRECISION X(N), F(N), D(N), WK(LWK) LOGICAL SPLINE C C DECLARE LOCAL VARIABLES. C INTEGER IC(2), INCFD DOUBLE PRECISION VC(2) DATA IC(1) /0/ DATA IC(2) /0/ DATA INCFD /1/ C C C***FIRST EXECUTABLE STATEMENT DPCHEZ C IF ( SPLINE ) THEN CALL DPCHSP (IC, VC, N, X, F, D, INCFD, WK, LWK, IERR) ELSE CALL DPCHIM (N, X, F, D, INCFD, IERR) ENDIF C C ERROR CONDITIONS ALREADY CHECKED IN DPCHSP OR DPCHIM RETURN C------------- LAST LINE OF DPCHEZ FOLLOWS ------------------------------ END SUBROUTINE DPCHIM(N,X,F,D,INCFD,IERR) C***BEGIN PROLOGUE DPCHIM C***DATE WRITTEN 811103 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E1B C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=DOUBLE PRECISION(DPCHIM-S DPCHIM-D), C CUBIC HERMITE INTERPOLATION,MONOTONE INTERPOLATION, C PIECEWISE CUBIC INTERPOLATION C***AUTHOR FRITSCH, F. N., (LLNL) C MATHEMATICS AND STATISTICS DIVISION C LAWRENCE LIVERMORE NATIONAL LABORATORY C P.O. BOX 808 (L-316) C LIVERMORE, CA 94550 C FTS 532-4275, (415) 422-4275 C***PURPOSE Set derivatives needed to determine a monotone piecewise C cubic Hermite interpolant to given data. Boundary values C are provided which are compatible with monotonicity. The C interpolant will have an extremum at each point where mono- C tonicity switches direction. (See DPCHIC if user control C is desired over boundary or switch conditions.) C***DESCRIPTION C C **** Double Precision version of DPCHIM **** C C DPCHIM: Piecewise Cubic Hermite Interpolation to C Monotone data. C C Sets derivatives needed to determine a monotone piecewise cubic C Hermite interpolant to the data given in X and F. C C Default boundary conditions are provided which are compatible C with monotonicity. (See DPCHIC if user control of boundary con- C ditions is desired.) C C If the data are only piecewise monotonic, the interpolant will C have an extremum at each point where monotonicity switches direc- C tion. (See DPCHIC if user control is desired in such cases.) C C To facilitate two-dimensional applications, includes an increment C between successive values of the F- and D-arrays. C C The resulting piecewise cubic Hermite function may be evaluated C by DPCHFE or DPCHFD. C C ---------------------------------------------------------------------- C C Calling sequence: C C PARAMETER (INCFD = ...) C INTEGER N, IERR C DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N) C C CALL DPCHIM (N, X, F, D, INCFD, IERR) C C Parameters: C C N -- (input) number of data points. (Error return if N.LT.2 .) C If N=2, simply does linear interpolation. C C X -- (input) real*8 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. C (Error return if not.) C C F -- (input) real*8 array of dependent variable values to be C interpolated. F(1+(I-1)*INCFD) is value corresponding to C X(I). DPCHIM is designed for monotonic data, but it will C work for any F-array. It will force extrema at points where C monotonicity switches direction. If some other treatment of C switch points is desired, DPCHIC should be used instead. C ----- C D -- (output) real*8 array of derivative values at the data C points. If the data are monotonic, these values will C determine a monotone cubic Hermite function. C The value corresponding to X(I) is stored in C D(1+(I-1)*INCFD), I=1(1)N. C No other entries in D are changed. C C INCFD -- (input) increment between successive values in F and D. C This argument is provided primarily for 2-D applications. C (Error return if INCFD.LT.1 .) C C IERR -- (output) error flag. C Normal return: C IERR = 0 (no errors). C Warning error: C IERR.GT.0 means that IERR switches in the direction C of monotonicity were detected. C "Recoverable" errors: C IERR = -1 if N.LT.2 . C IERR = -2 if INCFD.LT.1 . C IERR = -3 if the X-array is not strictly increasing. C (The D-array has not been changed in any of these cases.) C NOTE: The above errors are checked in the order listed, C and following arguments have **NOT** been validated. C C***REFERENCES 1. F.N.FRITSCH AND R.E.CARLSON, 'MONOTONE PIECEWISE C CUBIC INTERPOLATION,' SIAM J.NUMER.ANAL. 17, 2 (APRIL C 1980), 238-246. C 2. F.N.FRITSCH AND J.BUTLAND, 'A METHOD FOR CONSTRUCTING C LOCAL MONOTONE PIECEWISE CUBIC INTERPOLANTS,' LLNL C PREPRINT UCRL-87559 (APRIL 1982). C***ROUTINES CALLED DPCHST,XERROR C***END PROLOGUE DPCHIM C C ---------------------------------------------------------------------- C C Change record: C 82-02-01 1. Introduced DPCHST to reduce possible over/under- C flow problems. C 2. Rearranged derivative formula for same reason. C 82-06-02 1. Modified end conditions to be continuous functions C of data when monotonicity switches in next interval. C 2. Modified formulas so end conditions are less prone C of over/underflow problems. C 82-08-03 Minor cosmetic changes for release 1. C 87-07-07 Corrected XERROR calls for d.p. name(s). C C ---------------------------------------------------------------------- C C Programming notes: C C 1. The function DPCHST(ARG1,ARG2) is assumed to return zero if C either argument is zero, +1 if they are of the same sign, and C -1 if they are of opposite sign. C 2. To produce a single precision version, simply: C a. Change DPCHIM to PCHIM wherever it occurs, C b. Change DPCHST to PCHST wherever it occurs, C c. Change all references to the Fortran intrinsics to their C single precision equivalents, C d. Change the double precision declarations to real, and C e. Change the constants ZERO and THREE to single precision. C C DECLARE ARGUMENTS. C INTEGER N, INCFD, IERR DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N) C C DECLARE LOCAL VARIABLES. C INTEGER I, NLESS1 DOUBLE PRECISION DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE, * H1, H2, HSUM, HSUMT3, THREE, W1, W2, ZERO DOUBLE PRECISION DPCHST DATA ZERO /0.D0/, THREE/3.D0/ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT DPCHIM IF ( N.LT.2 ) GO TO 5001 IF ( INCFD.LT.1 ) GO TO 5002 DO 1 I = 2, N IF ( X(I).LE.X(I-1) ) GO TO 5003 1 CONTINUE C C FUNCTION DEFINITION IS OK, GO ON. C IERR = 0 NLESS1 = N - 1 H1 = X(2) - X(1) DEL1 = (F(1,2) - F(1,1))/H1 DSAVE = DEL1 C C SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. C IF (NLESS1 .GT. 1) GO TO 10 D(1,1) = DEL1 D(1,N) = DEL1 GO TO 5000 C C NORMAL CASE (N .GE. 3). C 10 CONTINUE H2 = X(3) - X(2) DEL2 = (F(1,3) - F(1,2))/H2 C C SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE C SHAPE-PRESERVING. C HSUM = H1 + H2 W1 = (H1 + HSUM)/HSUM W2 = -H1/HSUM D(1,1) = W1*DEL1 + W2*DEL2 IF ( DPCHST(D(1,1),DEL1) .LE. ZERO) THEN D(1,1) = ZERO ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO) THEN C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. DMAX = THREE*DEL1 IF (DABS(D(1,1)) .GT. DABS(DMAX)) D(1,1) = DMAX ENDIF C C LOOP THROUGH INTERIOR POINTS. C DO 50 I = 2, NLESS1 IF (I .EQ. 2) GO TO 40 C H1 = H2 H2 = X(I+1) - X(I) HSUM = H1 + H2 DEL1 = DEL2 DEL2 = (F(1,I+1) - F(1,I))/H2 40 CONTINUE C C SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. C D(1,I) = ZERO IF ( DPCHST(DEL1,DEL2) ) 42, 41, 45 C C COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY. C 41 CONTINUE IF (DEL2 .EQ. ZERO) GO TO 50 IF ( DPCHST(DSAVE,DEL2) .LT. ZERO) IERR = IERR + 1 DSAVE = DEL2 GO TO 50 C 42 CONTINUE IERR = IERR + 1 DSAVE = DEL2 GO TO 50 C C USE BRODLIE MODIFICATION OF BUTLAND FORMULA. C 45 CONTINUE HSUMT3 = HSUM+HSUM+HSUM W1 = (HSUM + H1)/HSUMT3 W2 = (HSUM + H2)/HSUMT3 DMAX = DMAX1( DABS(DEL1), DABS(DEL2) ) DMIN = DMIN1( DABS(DEL1), DABS(DEL2) ) DRAT1 = DEL1/DMAX DRAT2 = DEL2/DMAX D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2) C 50 CONTINUE C C SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE C SHAPE-PRESERVING. C W1 = -H2/HSUM W2 = (H2 + HSUM)/HSUM D(1,N) = W1*DEL1 + W2*DEL2 IF ( DPCHST(D(1,N),DEL2) .LE. ZERO) THEN D(1,N) = ZERO ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO) THEN C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. DMAX = THREE*DEL2 IF (DABS(D(1,N)) .GT. DABS(DMAX)) D(1,N) = DMAX ENDIF C C NORMAL RETURN. C 5000 CONTINUE RETURN C C ERROR RETURNS. C 5001 CONTINUE C N.LT.2 RETURN. IERR = -1 CALL XERROR ('DPCHIM -- NUMBER OF DATA POINTS LESS THAN TWO' * , 0, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('DPCHIM -- INCREMENT LESS THAN ONE' * , 0, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('DPCHIM -- X-ARRAY NOT STRICTLY INCREASING' * , 0, IERR, 1) RETURN C------------- LAST LINE OF DPCHIM FOLLOWS ----------------------------- END DOUBLE PRECISION FUNCTION DPCHST(ARG1,ARG2) C***BEGIN PROLOGUE DPCHST C***REFER TO DPCHCE,DPCHCI,DPCHCS,DPCHIM C***ROUTINES CALLED (NONE) C***DESCRIPTION C C DPCHST: DPCHIP Sign-Testing Routine. C C C Returns: C -1. if ARG1 and ARG2 are of opposite sign. C 0. if either argument is zero. C +1. if ARG1 and ARG2 are of the same sign. C C The object is to do this without multiplying ARG1*ARG2, to avoid C possible over/underflow problems. C C Fortran intrinsics used: SIGN. C C ---------------------------------------------------------------------- C C Programmed by: Fred N. Fritsch, FTS 532-4275, (415) 422-4275, C Mathematics and Statistics Division, C Lawrence Livermore National Laboratory. C C Change record: C 82-08-05 Converted to SLATEC library version. C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a single precision version, simply: C a. Change DPCHST to PCHST wherever it occurs, C b. Change all references to the Fortran intrinsics to their C single presision equivalents, C c. Change the double precision declarations to real, and C d. Change the constants ZERO and ONE to single precision. C***END PROLOGUE DPCHST C C DECLARE ARGUMENTS. C DOUBLE PRECISION ARG1, ARG2 C C DECLARE LOCAL VARIABLES. C DOUBLE PRECISION ONE, ZERO DATA ZERO /0.D0/, ONE/1.D0/ C C PERFORM THE TEST. C C***FIRST EXECUTABLE STATEMENT DPCHST DPCHST = DSIGN(ONE,ARG1) * DSIGN(ONE,ARG2) IF ((ARG1.EQ.ZERO) .OR. (ARG2.EQ.ZERO)) DPCHST = ZERO C RETURN C------------- LAST LINE OF DPCHST FOLLOWS ----------------------------- END SUBROUTINE DPCHSP(IC,VC,N,X,F,D,INCFD,WK,NWK,IERR) C***BEGIN PROLOGUE DPCHSP C***DATE WRITTEN 820503 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E1B C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=DOUBLE PRECISION(PCHSP-S DPCHSP-D), C CUBIC HERMITE INTERPOLATION,PIECEWISE CUBIC INTERPOLATION, C SPLINE INTERPOLATION C***AUTHOR FRITSCH, F. N., (LLNL) C MATHEMATICS AND STATISTICS DIVISION C LAWRENCE LIVERMORE NATIONAL LABORATORY C P.O. BOX 808 (L-316) C LIVERMORE, CA 94550 C FTS 532-4275, (415) 422-4275 C***PURPOSE Set derivatives needed to determine the Hermite represen- C tation of the cubic spline interpolant to given data, with C specified boundary conditions. C***DESCRIPTION C C **** Double Precision version of PCHSP **** C C DPCHSP: Piecewise Cubic Hermite Spline C C Computes the Hermite representation of the cubic spline inter- C polant to the data given in X and F satisfying the boundary C conditions specified by IC and VC. C C To facilitate two-dimensional applications, includes an increment C between successive values of the F- and D-arrays. C C The resulting piecewise cubic Hermite function may be evaluated C by DPCHFE or DPCHFD. C C NOTE: This is a modified version of C. de Boor'S cubic spline C routine CUBSPL. C C ---------------------------------------------------------------------- C C Calling sequence: C C PARAMETER (INCFD = ...) C INTEGER IC(2), N, NWK, IERR C DOUBLE PRECISION VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK) C C CALL DPCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR) C C Parameters: C C IC -- (input) integer array of length 2 specifying desired C boundary conditions: C IC(1) = IBEG, desired condition at beginning of data. C IC(2) = IEND, desired condition at end of data. C C IBEG = 0 to set D(1) so that the third derivative is con- C tinuous at X(2). This is the "not a knot" condition C provided by de Boor'S cubic spline routine CUBSPL. C < This is the default boundary condition. > C IBEG = 1 if first derivative at X(1) is given in VC(1). C IBEG = 2 if second derivative at X(1) is given in VC(1). C IBEG = 3 to use the 3-point difference formula for D(1). C (Reverts to the default b.c. if N.LT.3 .) C IBEG = 4 to use the 4-point difference formula for D(1). C (Reverts to the default b.c. if N.LT.4 .) C NOTES: C 1. An error return is taken if IBEG is out of range. C 2. For the "natural" boundary condition, use IBEG=2 and C VC(1)=0. C C IEND may take on the same values as IBEG, but applied to C derivative at X(N). In case IEND = 1 or 2, the value is C given in VC(2). C C NOTES: C 1. An error return is taken if IEND is out of range. C 2. For the "natural" boundary condition, use IEND=2 and C VC(2)=0. C C VC -- (input) real*8 array of length 2 specifying desired boundary C values, as indicated above. C VC(1) need be set only if IC(1) = 1 or 2 . C VC(2) need be set only if IC(2) = 1 or 2 . C C N -- (input) number of data points. (Error return if N.LT.2 .) C C X -- (input) real*8 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. C (Error return if not.) C C F -- (input) real*8 array of dependent variable values to be C interpolated. F(1+(I-1)*INCFD) is value corresponding to C X(I). C C D -- (output) real*8 array of derivative values at the data C points. These values will determine the cubic spline C interpolant with the requested boundary conditions. C The value corresponding to X(I) is stored in C D(1+(I-1)*INCFD), I=1(1)N. C No other entries in D are changed. C C INCFD -- (input) increment between successive values in F and D. C This argument is provided primarily for 2-D applications. C (Error return if INCFD.LT.1 .) C C WK -- (scratch) real*8 array of working storage. C C NWK -- (input) length of work array. C (Error return if NWK.LT.2*N .) C C IERR -- (output) error flag. C Normal return: C IERR = 0 (no errors). C "Recoverable" errors: C IERR = -1 if N.LT.2 . C IERR = -2 if INCFD.LT.1 . C IERR = -3 if the X-array is not strictly increasing. C IERR = -4 if IBEG.LT.0 or IBEG.GT.4 . C IERR = -5 if IEND.LT.0 of IEND.GT.4 . C IERR = -6 if both of the above are true. C IERR = -7 if NWK is too small. C NOTE: The above errors are checked in the order listed, C and following arguments have **NOT** been validated. C (The D-array has not been changed in any of these cases.) C IERR = -8 in case of trouble solving the linear system C for the interior derivative values. C (The D-array may have been changed in this case.) C ( Do **NOT** use it! ) C C***REFERENCES CARL DE BOOR, A PRACTICAL GUIDE TO SPLINES, SPRINGER- C VERLAG (NEW YORK, 1978), PP. 53-59. C***ROUTINES CALLED DPCHDF,XERROR C***END PROLOGUE DPCHSP C C ---------------------------------------------------------------------- C C Change record: C 82-08-04 Converted to SLATEC library version. C 87-07-07 Corrected XERROR calls for d.p. name(s). C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a single precision version, simply: C a. Change DPCHSP to PCHSP wherever it occurs, C b. Change the double precision declarations to real, and C c. Change the constants ZERO, HALF, ... to single precision. C C DECLARE ARGUMENTS. C INTEGER IC(2), N, INCFD, NWK, IERR DOUBLE PRECISION VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(2,N) C C DECLARE LOCAL VARIABLES. C INTEGER IBEG, IEND, INDEX, J, NM1 DOUBLE PRECISION G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), * ZERO DOUBLE PRECISION DPCHDF C DATA ZERO /0.D0/, HALF/.5D0/, ONE/1.D0/, TWO/2.D0/, THREE/3.D0/ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT DPCHSP IF ( N.LT.2 ) GO TO 5001 IF ( INCFD.LT.1 ) GO TO 5002 DO 1 J = 2, N IF ( X(J).LE.X(J-1) ) GO TO 5003 1 CONTINUE C IBEG = IC(1) IEND = IC(2) IERR = 0 IF ( (IBEG.LT.0).OR.(IBEG.GT.4) ) IERR = IERR - 1 IF ( (IEND.LT.0).OR.(IEND.GT.4) ) IERR = IERR - 2 IF ( IERR.LT.0 ) GO TO 5004 C C FUNCTION DEFINITION IS OK -- GO ON. C IF ( NWK .LT. 2*N ) GO TO 5007 C C COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO, C COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.). DO 5 J=2,N WK(1,J) = X(J) - X(J-1) WK(2,J) = (F(1,J) - F(1,J-1))/WK(1,J) 5 CONTINUE C C SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. C IF ( IBEG.GT.N ) IBEG = 0 IF ( IEND.GT.N ) IEND = 0 C C SET UP FOR BOUNDARY CONDITIONS. C IF ( (IBEG.EQ.1).OR.(IBEG.EQ.2) ) THEN D(1,1) = VC(1) ELSE IF (IBEG .GT. 2) THEN C PICK UP FIRST IBEG POINTS, IN REVERSE ORDER. DO 10 J = 1, IBEG INDEX = IBEG-J+1 C INDEX RUNS FROM IBEG DOWN TO 1. XTEMP(J) = X(INDEX) IF (J .LT. IBEG) STEMP(J) = WK(2,INDEX) 10 CONTINUE C -------------------------------- D(1,1) = DPCHDF (IBEG, XTEMP, STEMP, IERR) C -------------------------------- IF (IERR .NE. 0) GO TO 5009 IBEG = 1 ENDIF C IF ( (IEND.EQ.1).OR.(IEND.EQ.2) ) THEN D(1,N) = VC(2) ELSE IF (IEND .GT. 2) THEN C PICK UP LAST IEND POINTS. DO 15 J = 1, IEND INDEX = N-IEND+J C INDEX RUNS FROM N+1-IEND UP TO N. XTEMP(J) = X(INDEX) IF (J .LT. IEND) STEMP(J) = WK(2,INDEX+1) 15 CONTINUE C -------------------------------- D(1,N) = DPCHDF (IEND, XTEMP, STEMP, IERR) C -------------------------------- IF (IERR .NE. 0) GO TO 5009 IEND = 1 ENDIF C C --------------------( BEGIN CODING FROM CUBSPL )-------------------- C C **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF C F AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM- C INATION, WITH S(J) ENDING UP IN D(1,J), ALL J. C WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE. C C CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM C WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1) C IF (IBEG .EQ. 0) THEN IF (N .EQ. 2) THEN C NO CONDITION AT LEFT END AND N = 2. WK(2,1) = ONE WK(1,1) = ONE D(1,1) = TWO*WK(2,2) ELSE C 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 .EQ. 1) THEN C SLOPE PRESCRIBED AT LEFT END. WK(2,1) = ONE WK(1,1) = ZERO ELSE C 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 C C IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND C CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH C EQUATION READS WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J). C NM1 = N-1 IF (NM1 .GT. 1) THEN DO 20 J=2,NM1 IF (WK(2,J-1) .EQ. 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)) 20 CONTINUE ENDIF C C CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM C (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N) C C IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK- C SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT C AT THIS POINT. IF (IEND .EQ. 1) GO TO 30 C IF (IEND .EQ. 0) THEN IF (N.EQ.2 .AND. IBEG.EQ.0) THEN C 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.EQ.2) .OR. (N.EQ.3 .AND. IBEG.EQ.0)) THEN C EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT* C NOT-A-KNOT AT LEFT END POINT). D(1,N) = TWO*WK(2,N) WK(2,N) = ONE IF (WK(2,N-1) .EQ. ZERO) GO TO 5008 G = -ONE/WK(2,N-1) ELSE C NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR ALSO NOT-A- C KNOT AT LEFT END POINT. G = WK(1,N-1) + WK(1,N) C 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) .EQ. ZERO) GO TO 5008 G = -G/WK(2,N-1) WK(2,N) = WK(1,N-1) ENDIF ELSE C 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) .EQ. ZERO) GO TO 5008 G = -ONE/WK(2,N-1) ENDIF C C COMPLETE FORWARD PASS OF GAUSS ELIMINATION. C WK(2,N) = G*WK(1,N-1) + WK(2,N) IF (WK(2,N) .EQ. ZERO) GO TO 5008 D(1,N) = (G*D(1,N-1) + D(1,N))/WK(2,N) C C CARRY OUT BACK SUBSTITUTION C 30 CONTINUE DO 40 J=NM1,1,-1 IF (WK(2,J) .EQ. ZERO) GO TO 5008 D(1,J) = (D(1,J) - WK(1,J)*D(1,J+1))/WK(2,J) 40 CONTINUE C --------------------( END CODING FROM CUBSPL )-------------------- C C NORMAL RETURN. C RETURN C C ERROR RETURNS. C 5001 CONTINUE C N.LT.2 RETURN. IERR = -1 CALL XERROR ('DPCHSP -- NUMBER OF DATA POINTS LESS THAN TWO' * , 0, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('DPCHSP -- INCREMENT LESS THAN ONE' * , 0, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('DPCHSP -- X-ARRAY NOT STRICTLY INCREASING' * , 0, IERR, 1) RETURN C 5004 CONTINUE C IC OUT OF RANGE RETURN. IERR = IERR - 3 CALL XERROR ('DPCHSP -- IC OUT OF RANGE' * , 0, IERR, 1) RETURN C 5007 CONTINUE C NWK TOO SMALL RETURN. IERR = -7 CALL XERROR ('DPCHSP -- WORK ARRAY TOO SMALL' * , 0, IERR, 1) RETURN C 5008 CONTINUE C SINGULAR SYSTEM. C *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES *** C *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). *** IERR = -8 CALL XERROR ('DPCHSP -- SINGULAR LINEAR SYSTEM' * , 0, IERR, 1) RETURN C 5009 CONTINUE C ERROR RETURN FROM DPCHDF. C *** THIS CASE SHOULD NEVER OCCUR *** IERR = -9 CALL XERROR ('DPCHSP -- ERROR RETURN FROM DPCHDF' * , 0, IERR, 1) RETURN C------------- LAST LINE OF DPCHSP FOLLOWS ----------------------------- END DOUBLE PRECISION FUNCTION DPCHDF(K,X,S,IERR) C***BEGIN PROLOGUE DPCHDF C***REFER TO DPCHCE,DPCHSP C***ROUTINES CALLED XERROR C***REVISION DATE 870707 (YYMMDD) C***DESCRIPTION C C DPCHDF: DPCHIP Finite Difference Formula C C Uses a divided difference formulation to compute a K-point approx- C imation to the derivative at X(K) based on the data in X and S. C C Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary C derivative approximations. C C ---------------------------------------------------------------------- C C On input: C K is the order of the desired derivative approximation. C K must be at least 3 (error return if not). C X contains the K values of the independent variable. C X need not be ordered, but the values **MUST** be C distinct. (Not checked here.) C S contains the associated slope values: C S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. C (Note that S need only be of length K-1.) C C On return: C S will be destroyed. C IERR will be set to -1 if K.LT.2 . C DPCHDF will be set to the desired derivative approximation if C IERR=0 or to zero if IERR=-1. C C ---------------------------------------------------------------------- C C Reference: Carl de Boor, A Practical Guide to Splines, Springer- C Verlag (New York, 1978), pp. 10-16. C C***END PROLOGUE DPCHDF C C ---------------------------------------------------------------------- C C Programmed by: Fred N. Fritsch, FTS 532-4275, (415) 422-4275, C Mathematics and Statistics Division, C Lawrence Livermore National Laboratory. C C Change record: C 82-08-05 Converted to SLATEC library version. C 87-07-07 Corrected XERROR calls for d.p. name(s). C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a single precision version, simply: C a. Change DPCHDF to PCHDF wherever it occurs, C b. Change the double precision declarations to real, and C c. Change the constant Zero to single precision. C C DECLARE ARGUMENTS. C INTEGER K, IERR DOUBLE PRECISION X(K), S(K) C C DECLARE LOCAL VARIABLES. C INTEGER I, J DOUBLE PRECISION VALUE, ZERO DATA ZERO /0.D0/ C C CHECK FOR LEGAL VALUE OF K. C C***FIRST EXECUTABLE STATEMENT DPCHDF IF (K .LT. 3) GO TO 5001 C C COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. C DO 10 J = 2, K-1 DO 9 I = 1, K-J S(I) = (S(I+1)-S(I))/(X(I+J)-X(I)) 9 CONTINUE 10 CONTINUE C C EVALUATE DERIVATIVE AT X(K). C VALUE = S(1) DO 20 I = 2, K-1 VALUE = S(I) + VALUE*(X(K)-X(I)) 20 CONTINUE C C NORMAL RETURN. C IERR = 0 DPCHDF = VALUE RETURN C C ERROR RETURN. C 5001 CONTINUE C K.LT.3 RETURN. IERR = -1 CALL XERROR ('DPCHDF -- K LESS THAN THREE' * , 0, IERR, 1) DPCHDF = ZERO RETURN C------------- LAST LINE OF DPCHDF FOLLOWS ----------------------------- END SUBROUTINE DPCHEV(N,X,F,D,NVAL,XVAL,FVAL,DVAL,IERR) C***BEGIN PROLOGUE DPCHEV C***DATE WRITTEN 870828 (YYMMDD) C***REVISION DATE 870828 (YYMMDD) C***CATEGORY NO. E3,H1 C***KEYWORDS CUBIC HERMITE OR SPLINE DIFFERENTIATION,CUBIC HERMITE C EVALUATION,EASY TO USE SPLINE OR CUBIC HERMITE EVALUATOR C***AUTHOR KAHANER, D.K., (NBS) C SCIENTIFIC COMPUTING DIVISION C NATIONAL BUREAU OF STANDARDS C ROOM A161, TECHNOLOGY BUILDING C GAITHERSBURG, MARYLAND 20899 C (301) 975-3808 C***PURPOSE Evaluates the function and first derivative of a piecewise C cubic Hermite or spline function at an array of points C XVAL. It is easy to use. C***DESCRIPTION C C DPCHEV: Piecewise Cubic Hermite or Spline Derivative C Evaluator. Easy to Use. C C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C C Evaluates the function and first derivative of the cubic Hermite C or spline function defined by N, X, F, D, at the array of points C XVAL. C C C This is an easy to use driver for the routines by F.N. Fritsch C described in reference (2) below. Those also have other C capabilities. C C ---------------------------------------------------------------------- C C Calling sequence: CALL DPCHEV (N, X, F, D, NVAL, XVAL, FVAL, DVAL, IERR) C C INTEGER N, NVAL, IERR C DOUBLE PRECISION X(N), F(N), D(N), XVAL(NVAL), FVAL(NVAL), DVAL(NVAL) C C Parameters: C C N -- (input) number of data points. (Error return if N.LT.2 .) C C X -- (input) double precision array of independent variable C values. The elements of X must be strictly increasing: C X(I-1) .LT. X(I), I = 2(1)N. (Error return if not.) C C F -- (input) double precision array of function values. F(I) is C the value corresponding to X(I). C C D -- (input) double precision array of derivative values. C D(I) is the value corresponding to X(I). C C NVAL -- (input) number of points at which the functions are to be C evaluated. ( Error return if NVAL.LT.1 ) C C XVAL -- (input) double precision array of points at which the C functions are to be evaluated. C C NOTES: C 1. The evaluation will be most efficient if the elements C of XVAL are increasing relative to X; C that is, XVAL(J) .GE. X(I) C implies XVAL(K) .GE. X(I), all K.GE.J . C 2. If any of the XVAL are outside the interval [X(1),X(N)], C values are extrapolated from the nearest extreme cubic, C and a warning error is returned. C C FVAL -- (output) double precision array of values of the cubic C Hermite function defined by N, X, F, D at the points XVAL. C C DVAL -- (output) double precision array of values of the C first derivative of the same function at the points XVAL. C C IERR -- (output) error flag. C Normal return: C IERR = 0 (no errors). C Warning error: C IERR.GT.0 means that extrapolation was performed at C IERR points. C "Recoverable" errors: C IERR = -1 if N.LT.2 . C IERR = -3 if the X-array is not strictly increasing. C IERR = -4 if NVAL.LT.1 . C (Output arrays have not been changed in any of these cases.) C NOTE: The above errors are checked in the order listed, C and following arguments have **NOT** been validated. C IERR = -5 if an error has occurred in the lower-level C routine DCHFDV. NB: this should never happen. C Notify the author **IMMEDIATELY** if it does. C C ---------------------------------------------------------------------- C***REFERENCES 1. F.N.FRITSCH AND R.E.CARLSON, 'MONOTONE PIECEWISE C CUBIC INTERPOLATION,' SIAM J.NUMER.ANAL. 17, 2 (APRIL C 1980), 238-246. C 2. F.N.FRITSCH, 'PIECEWISE CUBIC HERMITE INTERPOLATION C PACKAGE, FINAL SPECIFICATIONS', LAWRENCE LIVERMORE C NATIONAL LABORATORY, COMPUTER DOCUMENTATION UCID-30194 C AUGUST 1982. C***ROUTINES CALLED DPCHFD C***END PROLOGUE DPCHEV INTEGER N, NVAL, IERR DOUBLE PRECISION X(N), F(N), D(N), XVAL(NVAL), FVAL(NVAL), *DVAL(NVAL) C C DECLARE LOCAL VARIABLES. C INTEGER INCFD LOGICAL SKIP DATA SKIP /.TRUE./ DATA INCFD /1/ C C C***FIRST EXECUTABLE STATEMENT DPCHEV C CALL DPCHFD(N,X,F,D,INCFD,SKIP,NVAL,XVAL,FVAL,DVAL,IERR) C C 5000 CONTINUE RETURN C C------------- LAST LINE OF DPCHEV FOLLOWS ------------------------------ END SUBROUTINE DPCHFD(N,X,F,D,INCFD,SKIP,NE,XE,FE,DE,IERR) C***BEGIN PROLOGUE DPCHFD C***DATE WRITTEN 811020 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E3,H1 C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=DOUBLE PRECISION(PCHFD-S DPCHFD-D), C CUBIC HERMITE DIFFERENTIATION,CUBIC HERMITE EVALUATION, C HERMITE INTERPOLATION,PIECEWISE CUBIC EVALUATION C***AUTHOR FRITSCH, F. N., (LLNL) C MATHEMATICS AND STATISTICS DIVISION C LAWRENCE LIVERMORE NATIONAL LABORATORY C P.O. BOX 808 (L-316) C LIVERMORE, CA 94550 C FTS 532-4275, (415) 422-4275 C***PURPOSE Evaluate a piecewise cubic hermite function and its first C derivative at an array of points. May be used by itself C for Hermite interpolation, or as an evaluator for DPCHIM C or DPCHIC. If only function values are required, use C DPCHFE instead. C***DESCRIPTION C C **** Double Precision version of PCHFD **** C C DPCHFD: Piecewise Cubic Hermite Function and Derivative C evaluator C C Evaluates the cubic Hermite function defined by N, X, F, D, to- C gether with its first derivative, at the points XE(J), J=1(1)NE. C C If only function values are required, use DPCHFE, instead. C C To provide compatibility with DPCHIM and DPCHIC, includes an C increment between successive values of the F- and D-arrays. C C ---------------------------------------------------------------------- C C Calling sequence: C C PARAMETER (INCFD = ...) C INTEGER N, NE, IERR C DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE), C DE(NE) C LOGICAL SKIP C C CALL DPCHFD (N, X, F, D, INCFD, SKIP, NE, XE, FE, DE, IERR) C C Parameters: C C N -- (input) number of data points. (Error return if N.LT.2 .) C C X -- (input) real*8 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. C (Error return if not.) C C F -- (input) real*8 array of function values. F(1+(I-1)*INCFD) is C the value corresponding to X(I). C C D -- (input) real*8 array of derivative values. D(1+(I-1)*INCFD) C is the value corresponding to X(I). C C INCFD -- (input) increment between successive values in F and D. C (Error return if INCFD.LT.1 .) C C SKIP -- (input/output) logical variable which should be set to C .TRUE. if the user wishes to skip checks for validity of C preceding parameters, or to .FALSE. otherwise. C This will save time in case these checks have already C been performed (say, in DPCHIM or DPCHIC). C SKIP will be set to .TRUE. on normal return. C C NE -- (input) number of evaluation points. (Error return if C NE.LT.1 .) C C XE -- (input) real*8 array of points at which the functions are to C be evaluated. C C C NOTES: C 1. The evaluation will be most efficient if the elements C of XE are increasing relative to X; C that is, XE(J) .GE. X(I) C implies XE(K) .GE. X(I), all K.GE.J . C 2. If any of the XE are outside the interval [X(1),X(N)], C values are extrapolated from the nearest extreme cubic, C and a warning error is returned. C C FE -- (output) real*8 array of values of the cubic Hermite C function defined by N, X, F, D at the points XE. C C DE -- (output) real*8 array of values of the first derivative of C the same function at the points XE. C C IERR -- (output) error flag. C Normal return: C IERR = 0 (no errors). C Warning error: C IERR.GT.0 means that extrapolation was performed at C IERR points. C "Recoverable" errors: C IERR = -1 if N.LT.2 . C IERR = -2 if INCFD.LT.1 . C IERR = -3 if the X-array is not strictly increasing. C IERR = -4 if NE.LT.1 . C (Output arrays have not been changed in any of these cases.) C NOTE: The above errors are checked in the order listed, C and following arguments have **NOT** been validated. C IERR = -5 if an error has occurred in the lower-level C routine DCHFDV. NB: this should never happen. C Notify the author **IMMEDIATELY** if it does. C C***REFERENCES (NONE) C***ROUTINES CALLED DCHFDV,XERROR C***END PROLOGUE DPCHFD C C ---------------------------------------------------------------------- C C Change record: C 82-08-03 Minor cosmetic changes for release 1. C 87-07-07 Corrected XERROR calls for d.p. name(s). C C ---------------------------------------------------------------------- C C Programming notes: C C 1. To produce a single precision version, simply: C a. Change DPCHFD to PCHFD, and DCHFDV to CHFDV, wherever they C occur, C b. Change the double precision declaration to real, C C 2. Most of the coding between the call to DCHFDV and the end of C the IR-loop could be eliminated if it were permissible to C assume that XE is ordered relative to X. C C 3. DCHFDV does not assume that X1 is less than X2. thus, it would C be possible to write a version of DPCHFD that assumes a strict- C ly decreasing X-array by simply running the IR-loop backwards C (and reversing the order of appropriate tests). C C 4. The present code has a minor bug, which I have decided is not C worth the effort that would be required to fix it. C If XE contains points in [X(N-1),X(N)], followed by points .LT. C X(N-1), followed by points .GT.X(N), the extrapolation points C will be counted (at least) twice in the total returned in IERR. C C DECLARE ARGUMENTS. C INTEGER N, INCFD, NE, IERR DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE), * DE(NE) LOGICAL SKIP C C DECLARE LOCAL VARIABLES. C INTEGER I, IERC, IR, J, JFIRST, NEXT(2), NJ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT DPCHFD IF (SKIP) GO TO 5 C IF ( N.LT.2 ) GO TO 5001 IF ( INCFD.LT.1 ) GO TO 5002 DO 1 I = 2, N IF ( X(I).LE.X(I-1) ) GO TO 5003 1 CONTINUE C C FUNCTION DEFINITION IS OK, GO ON. C 5 CONTINUE IF ( NE.LT.1 ) GO TO 5004 IERR = 0 SKIP = .TRUE. C C LOOP OVER INTERVALS. ( INTERVAL INDEX IS IL = IR-1 . ) C ( INTERVAL IS X(IL).LE.X.LT.X(IR) . ) JFIRST = 1 IR = 2 10 CONTINUE C C SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS. C IF (JFIRST .GT. NE) GO TO 5000 C C LOCATE ALL POINTS IN INTERVAL. C DO 20 J = JFIRST, NE IF (XE(J) .GE. X(IR)) GO TO 30 20 CONTINUE J = NE + 1 GO TO 40 C C HAVE LOCATED FIRST POINT BEYOND INTERVAL. C 30 CONTINUE IF (IR .EQ. N) J = NE + 1 C 40 CONTINUE NJ = J - JFIRST C C SKIP EVALUATION IF NO POINTS IN INTERVAL. C IF (NJ .EQ. 0) GO TO 50 C C EVALUATE CUBIC AT XE(I), I = JFIRST (1) J-1 . C C ---------------------------------------------------------------- CALL DCHFDV (X(IR-1),X(IR), F(1,IR-1),F(1,IR), D(1,IR-1),D(1,IR) * ,NJ, XE(JFIRST), FE(JFIRST), DE(JFIRST), NEXT, IERC) C ---------------------------------------------------------------- IF (IERC .LT. 0) GO TO 5005 C IF (NEXT(2) .EQ. 0) GO TO 42 C IF (NEXT(2) .GT. 0) THEN C IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE C RIGHT OF X(IR). C IF (IR .LT. N) GO TO 41 C IF (IR .EQ. N) THEN C THESE ARE ACTUALLY EXTRAPOLATION POINTS. IERR = IERR + NEXT(2) GO TO 42 41 CONTINUE C ELSE C WE SHOULD NEVER HAVE GOTTEN HERE. GO TO 5005 C ENDIF C ENDIF 42 CONTINUE C IF (NEXT(1) .EQ. 0) GO TO 49 C IF (NEXT(1) .GT. 0) THEN C IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE C LEFT OF X(IR-1). C IF (IR .GT. 2) GO TO 43 C IF (IR .EQ. 2) THEN C THESE ARE ACTUALLY EXTRAPOLATION POINTS. IERR = IERR + NEXT(1) GO TO 49 43 CONTINUE C ELSE C XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST C EVALUATION INTERVAL. C C FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1). DO 44 I = JFIRST, J-1 IF (XE(I) .LT. X(IR-1)) GO TO 45 44 CONTINUE C NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR C IN DCHFDV. GO TO 5005 C 45 CONTINUE C RESET J. (THIS WILL BE THE NEW JFIRST.) J = I C C NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY. DO 46 I = 1, IR-1 IF (XE(J) .LT. X(I)) GO TO 47 46 CONTINUE C NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1). C 47 CONTINUE C AT THIS POINT, EITHER XE(J) .LT. X(1) C OR X(I-1) .LE. XE(J) .LT. X(I) . C RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE C CYCLING. IR = MAX0(1, I-1) C ENDIF C ENDIF 49 CONTINUE C JFIRST = J C C END OF IR-LOOP. C 50 CONTINUE IR = IR + 1 IF (IR .LE. N) GO TO 10 C C NORMAL RETURN. C 5000 CONTINUE RETURN C C ERROR RETURNS. C 5001 CONTINUE C N.LT.2 RETURN. IERR = -1 CALL XERROR ('DPCHFD -- NUMBER OF DATA POINTS LESS THAN TWO' * , 0, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('DPCHFD -- INCREMENT LESS THAN ONE' * , 0, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('DPCHFD -- X-ARRAY NOT STRICTLY INCREASING' * , 0, IERR, 1) RETURN C 5004 CONTINUE C NE.LT.1 RETURN. IERR = -4 CALL XERROR ('DPCHFD -- NUMBER OF EVALUATION POINTS LESS THAN ONE' * , 0, IERR, 1) RETURN C 5005 CONTINUE C ERROR RETURN FROM DCHFDV. C *** THIS CASE SHOULD NEVER OCCUR *** IERR = -5 CALL XERROR ('DPCHFD -- ERROR RETURN FROM DCHFDV -- FATAL' * , 0, IERR, 2) RETURN C------------- LAST LINE OF DPCHFD FOLLOWS ----------------------------- END SUBROUTINE DCHFDV(X1,X2,F1,F2,D1,D2,NE,XE,FE,DE,NEXT,IERR) C***BEGIN PROLOGUE DCHFDV C***DATE WRITTEN 811019 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E3,H1 C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=DOUBLE PRECISION(CHFDV-S DCHFDV-D), C CUBIC HERMITE DIFFERENTIATION,CUBIC HERMITE EVALUATION, C CUBIC POLYNOMIAL EVALUATION C***AUTHOR FRITSCH, F. N., (LLNL) C MATHEMATICS AND STATISTICS DIVISION C LAWRENCE LIVERMORE NATIONAL LABORATORY C P.O. BOX 808 (L-316) C LIVERMORE, CA 94550 C FTS 532-4275, (415) 422-4275 C***PURPOSE Evaluate a cubic polynomial given in Hermite form and its C first derivative at an array of points. While designed for C use by DPCHFD, it may be useful directly as an evaluator C for a piecewise cubic Hermite function in applications, C such as graphing, where the interval is known in advance. C If only function values are required, use DCHFEV instead. C***DESCRIPTION C C **** Double Precision version of CHFDV **** C C DCHFDV: Cubic Hermite Function and Derivative Evaluator C C Evaluates the cubic polynomial determined by function values C F1,F2 and derivatives D1,D2 on interval (X1,X2), together with C its first derivative, at the points XE(J), J=1(1)NE. C C If only function values are required, use DCHFEV, instead. C C ---------------------------------------------------------------------- C C Calling sequence: C C INTEGER NE, NEXT(2), IERR C DOUBLE PRECISION X1, X2, F1, F2, D1, D2, XE(NE), FE(NE), C DE(NE) C C CALL DCHFDV (X1,X2, F1,F2, D1,D2, NE, XE, FE, DE, NEXT, IERR) C C Parameters: C C X1,X2 -- (input) endpoints of interval of definition of cubic. C (Error return if X1.EQ.X2 .) C C F1,F2 -- (input) values of function at X1 and X2, respectively. C C D1,D2 -- (input) values of derivative at X1 and X2, respectively. C C NE -- (input) number of evaluation points. (Error return if C NE.LT.1 .) C C XE -- (input) real*8 array of points at which the functions are to C be evaluated. If any of the XE are outside the interval C [X1,X2], a warning error is returned in NEXT. C C FE -- (output) real*8 array of values of the cubic function C defined by X1,X2, F1,F2, D1,D2 at the points XE. C C DE -- (output) real*8 array of values of the first derivative of C the same function at the points XE. C C NEXT -- (output) integer array indicating number of extrapolation C points: C NEXT(1) = number of evaluation points to left of interval. C NEXT(2) = number of evaluation points to right of interval. C C IERR -- (output) error flag. C Normal return: C IERR = 0 (no errors). C "Recoverable" errors: C IERR = -1 if NE.LT.1 . C IERR = -2 if X1.EQ.X2 . C (Output arrays have not been changed in either case.) C C***REFERENCES (NONE) C***ROUTINES CALLED XERROR C***END PROLOGUE DCHFDV C C ---------------------------------------------------------------------- C C Change record: C 82-08-03 Minor cosmetic changes for release 1. C 87-07-07 Corrected XERROR calls for d.p. names(s). C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a single precision version, simply: C a. Change DCHFDV to CHFDV wherever it occurs, C b. Change the double precision declaration to real, C c. Change the constant Zero to single precision, and C d. Change the names of the Fortran functions: AMAX1, AMIN1. C C DECLARE ARGUMENTS. C INTEGER NE, NEXT(2), IERR DOUBLE PRECISION X1, X2, F1, F2, D1, D2, XE(NE), FE(NE), DE(NE) C C DECLARE LOCAL VARIABLES. C INTEGER I DOUBLE PRECISION C2, C2T2, C3, C3T3, DEL1, DEL2, DELTA, H, X, * XMI, XMA, ZERO DATA ZERO /0.D0/ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT DCHFDV IF (NE .LT. 1) GO TO 5001 H = X2 - X1 IF (H .EQ. ZERO) GO TO 5002 C C INITIALIZE. C IERR = 0 NEXT(1) = 0 NEXT(2) = 0 XMI = DMIN1(ZERO, H) XMA = DMAX1(ZERO, H) C C COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1). C DELTA = (F2 - F1)/H DEL1 = (D1 - DELTA)/H DEL2 = (D2 - DELTA)/H C (DELTA IS NO LONGER NEEDED.) C2 = -(DEL1+DEL1 + DEL2) C2T2 = C2 + C2 C3 = (DEL1 + DEL2)/H C (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.) C3T3 = C3+C3+C3 C C EVALUATION LOOP. C DO 500 I = 1, NE X = XE(I) - X1 FE(I) = F1 + X*(D1 + X*(C2 + X*C3)) DE(I) = D1 + X*(C2T2 + X*C3T3) C COUNT EXTRAPOLATION POINTS. IF ( X.LT.XMI ) NEXT(1) = NEXT(1) + 1 IF ( X.GT.XMA ) NEXT(2) = NEXT(2) + 1 C (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.) 500 CONTINUE C C NORMAL RETURN. C RETURN C C ERROR RETURNS. C 5001 CONTINUE C NE.LT.1 RETURN. IERR = -1 CALL XERROR ('DCHFDV -- NUMBER OF EVALUATION POINTS LESS THAN ONE' * , 51, IERR, 1) RETURN C 5002 CONTINUE C X1.EQ.X2 RETURN. IERR = -2 CALL XERROR ('DCHFDV -- INTERVAL ENDPOINTS EQUAL' * , 34, IERR, 1) RETURN C------------- LAST LINE OF DCHFDV FOLLOWS ----------------------------- END DOUBLE PRECISION FUNCTION DPCHQA(N,X,F,D,A,B,IERR) C***BEGIN PROLOGUE DPCHQA C***DATE WRITTEN 870829 (YYMMDD) C***REVISION DATE 870829 (YYMMDD) C***CATEGORY NO. E3,H2A2 C***KEYWORDS EASY TO USE CUBIC HERMITE OR SPLINE INTEGRATION C NUMERICAL INTEGRATION, QUADRATURE C***AUTHOR KAHANER, D.K., (NBS) C SCIENTIFIC COMPUTING DIVISION C NATIONAL BUREAU OF STANDARDS C ROOM A161, TECHNOLOGY BUILDING C GAITHERSBURG, MARYLAND 20899 C (301) 975-3808 C***PURPOSE Evaluates the definite integral of a piecewise cubic Hermit C or spline function over an arbitrary interval, easy to use. C***DESCRIPTION C C DPCHQA: Piecewise Cubic Hermite or Spline Integrator, C Arbitrary limits, Easy to Use. C C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C C Evaluates the definite integral of the cubic Hermite or spline C function defined by N, X, F, D over the interval [A, B]. This C is an easy to use driver for the routine DPCHIA by F.N. Fritsch C described in reference (2) below. That routine also has other C capabilities. C ---------------------------------------------------------------------- C C Calling sequence: C C VALUE = DPCHQA (N, X, F, D, A, B, IERR) C C INTEGER N, IERR C DOUBLE PRECISION X(N), F(N), D(N), A, B C C Parameters: C C VALUE -- (output) VALUE of the requested integral. C C N -- (input) number of data points. (Error return if N.LT.2 .) C C X -- (input) double precision array of independent variable C values. The elements of X must be strictly increasing: C X(I-1) .LT. X(I), I = 2(1)N. C (Error return if not.) C C F -- (input) double precision array of function values. C F(I) is the value corresponding to X(I). C C D -- (input) double precision array of derivative values. D(I) is C the value corresponding to X(I). C C A,B -- (input) the limits of integration. C NOTE: There is no requirement that [A,B] be contained in C [X(1),X(N)]. However, the resulting integral value C will be highly suspect, if not. C C IERR -- (output) error flag. C Normal return: C IERR = 0 (no errors). C Warning errors: C IERR = 1 if A is outside the interval [X(1),X(N)]. C IERR = 2 if B is outside the interval [X(1),X(N)]. C IERR = 3 if both of the above are true. (Note that this C means that either [A,B] contains data interval C or the intervals do not intersect at all.) C "Recoverable" errors: C IERR = -1 if N.LT.2 . C IERR = -3 if the X-array is not strictly increasing. C (Value has not been computed in any of these cases.) C NOTE: The above errors are checked in the order listed, C and following arguments have **NOT** been validated. C C***REFERENCES 1. F.N.FRITSCH AND R.E.CARLSON, 'MONOTONE PIECEWISE C CUBIC INTERPOLATION,' SIAM J.NUMER.ANAL. 17, 2 (APRIL C 1980), 238-246. C 2. F.N.FRITSCH, 'PIECEWISE CUBIC HERMITE INTERPOLATION C PACKAGE, FINAL SPECIFICATIONS', LAWRENCE LIVERMORE C NATIONAL LABORATORY, COMPUTER DOCUMENTATION UCID-30194 C AUGUST 1982. C***ROUTINES CALLED DPCHIA C***END PROLOGUE DPCHQA INTEGER N, IERR DOUBLE PRECISION X(N), F(N), D(N), A, B C C DECLARE LOCAL VARIABLES. C INTEGER INCFD DOUBLE PRECISION DPCHIA LOGICAL SKIP C C INITIALIZE. C DATA INCFD /1/ DATA SKIP /.TRUE./ C C C***FIRST EXECUTABLE STATEMENT DPCHQA DPCHQA = DPCHIA( N, X, F, D, INCFD, SKIP, A, B, IERR ) C C ERROR MESSAGES ARE FROM LOWER LEVEL ROUTINES RETURN C C------------- LAST LINE OF DPCHQA FOLLOWS ------------------------------ END DOUBLE PRECISION FUNCTION DPCHIA(N,X,F,D,INCFD,SKIP,A,B,IERR) C***BEGIN PROLOGUE DPCHIA C***DATE WRITTEN 820730 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E3,H2A2 C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=DOUBLE PRECISION(PCHIA-S DPCHIA-D), C CUBIC HERMITE INTERPOLATION,NUMERICAL INTEGRATION, C QUADRATURE C***AUTHOR FRITSCH, F. N., (LLNL) C MATHEMATICS AND STATISTICS DIVISION C LAWRENCE LIVERMORE NATIONAL LABORATORY C P.O. BOX 808 (L-316) C LIVERMORE, CA 94550 C FTS 532-4275, (415) 422-4275 C***PURPOSE Evaluate the definite integral of a piecewise cubic C Hermite function over an arbitrary interval. C***DESCRIPTION C C **** Double Precision version of PCHIA **** C C DPCHIA: Piecewise Cubic Hermite Integrator, Arbitrary limits C C Evaluates the definite integral of the cubic Hermite function C defined by N, X, F, D over the interval [A, B]. C C To provide compatibility with DPCHIM and DPCHIC, includes an C increment between successive values of the F- and D-arrays. C C ---------------------------------------------------------------------- C C Calling sequence: C C PARAMETER (INCFD = ...) C INTEGER N, IERR C DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N), A, B C DOUBLE PRECISION VALUE, DPCHIA C LOGICAL SKIP C C VALUE = DPCHIA (N, X, F, D, INCFD, SKIP, A, B, IERR) C C Parameters: C C VALUE -- (output) VALUE of the requested integral. C C N -- (input) number of data points. (Error return if N.LT.2 .) C C X -- (input) real*8 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. C (Error return if not.) C C F -- (input) real*8 array of function values. F(1+(I-1)*INCFD) is C the value corresponding to X(I). C C D -- (input) real*8 array of derivative values. D(1+(I-1)*INCFD) C is the value corresponding to X(I). C C INCFD -- (input) increment between successive values in F and D. C (Error return if INCFD.LT.1 .) C C SKIP -- (input/output) logical variable which should be set to C .TRUE. if the user wishes to skip checks for validity of C preceding parameters, or to .FALSE. otherwise. C This will save time in case these checks have already C been performed (say, in DPCHIM or DPCHIC). C SKIP will be set to .TRUE. on return with IERR.GE.0 . C C A,B -- (input) the limits of integration. C NOTE: There is no requirement that [A,B] be contained in C [X(1),X(N)]. However, the resulting integral value C will be highly suspect, if not. C C IERR -- (output) error flag. C Normal return: C IERR = 0 (no errors). C Warning errors: C IERR = 1 if A is outside the interval [X(1),X(N)]. C IERR = 2 if B is outside the interval [X(1),X(N)]. C IERR = 3 if both of the above are true. (Note that this C means that either [A,B] contains data interval C or the intervals do not intersect at all.) C "Recoverable" errors: C IERR = -1 if N.LT.2 . C IERR = -2 if INCFD.LT.1 . C IERR = -3 if the X-array is not strictly increasing. C (Value has not been computed in any of these cases.) C NOTE: The above errors are checked in the order listed, C and following arguments have **NOT** been validated. C C***REFERENCES (NONE) C***ROUTINES CALLED DCHFIV,DPCHID,XERROR C***END PROLOGUE DPCHIA C C ---------------------------------------------------------------------- C C Change record: C 82-08-04 Converted to SLATEC library version. C 87-07-07 Corrected conversion to double precision. C 87-07-07 Corrected XERROR calls for d.p. name(s). C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a single precision version, simply: C a. Change DPCHIA to PCHIA wherever it occurs, C b. Change DPCHID to PCHID wherever it occurs, C c. Change DPCHIV to PCHIV wherever it occurs, C d. Change the double precision declarations to real, and C e. Change the constant ZERO to single precision. C C DECLARE ARGUMENTS. C INTEGER N, INCFD, IERR DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N), A, B LOGICAL SKIP C C DECLARE LOCAL VARIABLES. C INTEGER I, IA, IB, IERD, IERV, IL, IR DOUBLE PRECISION VALUE, XA, XB, ZERO DOUBLE PRECISION DCHFIV, DPCHID C C INITIALIZE. C DATA ZERO /0.D0/ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT DPCHIA IF (SKIP) GO TO 5 C IF ( N.LT.2 ) GO TO 5001 IF ( INCFD.LT.1 ) GO TO 5002 DO 1 I = 2, N IF ( X(I).LE.X(I-1) ) GO TO 5003 1 CONTINUE C C FUNCTION DEFINITION IS OK, GO ON. C 5 CONTINUE SKIP = .TRUE. IERR = 0 IF ( (A.LT.X(1)) .OR. (A.GT.X(N)) ) IERR = IERR + 1 IF ( (B.LT.X(1)) .OR. (B.GT.X(N)) ) IERR = IERR + 2 C C COMPUTE INTEGRAL VALUE. C IF (A .EQ. B) THEN VALUE = ZERO ELSE XA = DMIN1 (A, B) XB = DMAX1 (A, B) IF (XB .LE. X(2)) THEN C INTERVAL IS TO LEFT OF X(2), SO USE FIRST CUBIC. C -------------------------------------------- VALUE = DCHFIV (X(1),X(2), F(1,1),F(1,2), * D(1,1),D(1,2), A, B, IERV) C -------------------------------------------- IF (IERV .LT. 0) GO TO 5004 ELSE IF (XA .GE. X(N-1)) THEN C INTERVAL IS TO RIGHT OF X(N-1), SO USE LAST CUBIC. C ----------------------------------------------- VALUE = DCHFIV(X(N-1),X(N), F(1,N-1),F(1,N), * D(1,N-1),D(1,N), A, B, IERV) C ----------------------------------------------- IF (IERV .LT. 0) GO TO 5004 ELSE C 'NORMAL' CASE -- XA.LT.XB, XA.LT.X(N-1), XB.GT.X(2). C ......LOCATE IA AND IB SUCH THAT C X(IA-1).LT.XA.LE.X(IA).LE.X(IB).LE.XB.LE.X(IB+1) IA = 1 DO 10 I = 1, N-1 IF (XA .GT. X(I)) IA = I + 1 10 CONTINUE C IA = 1 IMPLIES XA.LT.X(1) . OTHERWISE, C IA IS LARGEST INDEX SUCH THAT X(IA-1).LT.XA,. C IB = N DO 20 I = N, IA, -1 IF (XB .LT. X(I)) IB = I - 1 20 CONTINUE C IB = N IMPLIES XB.GT.X(N) . OTHERWISE, C IB IS SMALLEST INDEX SUCH THAT XB.LT.X(IB+1) . C C ......COMPUTE THE INTEGRAL. IERV = 0 IF (IB .LT. IA) THEN C THIS MEANS IB = IA-1 AND C (A,B) IS A SUBSET OF (X(IB),X(IA)). C ------------------------------------------------ VALUE = DCHFIV (X(IB),X(IA), F(1,IB),F(1,IA), * D(1,IB),D(1,IA), A, B, IERV) C ------------------------------------------------ IF (IERV .LT. 0) GO TO 5004 ELSE C C FIRST COMPUTE INTEGRAL OVER (X(IA),X(IB)). IF (IB .EQ. IA) THEN VALUE = ZERO ELSE C --------------------------------------------- VALUE = DPCHID (N, X, F, D, INCFD, SKIP, IA, IB, IERD) C --------------------------------------------- IF (IERD .LT. 0) GO TO 5005 ENDIF C C THEN ADD ON INTEGRAL OVER (XA,X(IA)). IF (XA .LT. X(IA)) THEN IL = MAX0 (1, IA-1) IR = IL + 1 C ------------------------------------- VALUE = VALUE + DCHFIV (X(IL),X(IR), F(1,IL),F(1,IR), * D(1,IL),D(1,IR), XA, X(IA), IERV) C ------------------------------------- IF (IERV .LT. 0) GO TO 5004 ENDIF C C THEN ADD ON INTEGRAL OVER (X(IB),XB). IF (XB .GT. X(IB)) THEN IR = MIN0 (IB+1, N) IL = IR - 1 C ------------------------------------- VALUE = VALUE + DCHFIV (X(IL),X(IR), F(1,IL),F(1,IR), * D(1,IL),D(1,IR), X(IB), XB, IERV) C ------------------------------------- IF (IERV .LT. 0) GO TO 5004 ENDIF C C FINALLY, ADJUST SIGN IF NECESSARY. IF (A .GT. B) VALUE = -VALUE ENDIF ENDIF ENDIF C C NORMAL RETURN. C DPCHIA = VALUE RETURN C C ERROR RETURNS. C 5001 CONTINUE C N.LT.2 RETURN. IERR = -1 CALL XERROR ('DPCHIA -- NUMBER OF DATA POINTS LESS THAN TWO' * , 0, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('DPCHIA -- INCREMENT LESS THAN ONE' * , 0, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('DPCHIA -- X-ARRAY NOT STRICTLY INCREASING' * , 0, IERR, 1) RETURN C 5004 CONTINUE C TROUBLE IN DCHFIV. (SHOULD NEVER OCCUR.) IERR = -4 CALL XERROR ('DPCHIA -- TROUBLE IN DCHFIV' * , 0, IERR, 1) RETURN C 5005 CONTINUE C TROUBLE IN DPCHID. (SHOULD NEVER OCCUR.) IERR = -5 CALL XERROR ('DPCHIA -- TROUBLE IN DPCHID' * , 0, IERR, 1) RETURN C------------- LAST LINE OF DPCHIA FOLLOWS ----------------------------- END DOUBLE PRECISION FUNCTION DPCHID(N,X,F,D,INCFD,SKIP,IA,IB,IERR) C***BEGIN PROLOGUE DPCHID C***DATE WRITTEN 820723 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E1B,H2A2 C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=DOUBLE PRECISION(PCHID-S DPCHID-D), C CUBIC HERMITE INTERPOLATION,NUMERICAL INTEGRATION, C QUADRATURE C***AUTHOR FRITSCH, F. N., (LLNL) C MATHEMATICS AND STATISTICS DIVISION C LAWRENCE LIVERMORE NATIONAL LABORATORY C P.O. BOX 808 (L-316) C LIVERMORE, CA 94550 C FTS 532-4275, (415) 422-4275 C***PURPOSE Evaluate the definite integral of a piecewise cubic C Hermite function over an interval whose endpoints are data C points. C***DESCRIPTION C C **** Double Precision version of PCHID **** C C DPCHID: Piecewise Cubic Hermite Integrator, Data limits C C Evaluates the definite integral of the cubic Hermite function C defined by N, X, F, D over the interval [X(IA), X(IB)]. C C To provide compatibility with DPCHIM and DPCHIC, includes an C increment between successive values of the F- and D-arrays. C C ---------------------------------------------------------------------- C C Calling sequence: C C PARAMETER (INCFD = ...) C INTEGER N, IA, IB, IERR C DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N) C LOGICAL SKIP C C VALUE = DPCHID (N, X, F, D, INCFD, SKIP, IA, IB, IERR) C C Parameters: C C VALUE -- (output) VALUE of the requested integral. C C N -- (input) number of data points. (Error return if N.LT.2 .) C C X -- (input) real*8 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. C (Error return if not.) C C F -- (input) real*8 array of function values. F(1+(I-1)*INCFD) is C the value corresponding to X(I). C C D -- (input) real*8 array of derivative values. D(1+(I-1)*INCFD) C is the value corresponding to X(I). C C INCFD -- (input) increment between successive values in F and D. C (Error return if INCFD.LT.1 .) C C SKIP -- (input/output) logical variable which should be set to C .TRUE. if the user wishes to skip checks for validity of C preceding parameters, or to .FALSE. otherwise. C This will save time in case these checks have already C been performed (say, in DPCHIM or DPCHIC). C SKIP will be set to .TRUE. on return with IERR = 0 or -4. C C IA,IB -- (input) indices in X-array for the limits of integration. C both must be in the range [1,N]. (Error return if not.) C No restrictions on their relative values. C C IERR -- (output) error flag. C Normal return: C IERR = 0 (no errors). C "Recoverable" errors: C IERR = -1 if N.LT.2 . C IERR = -2 if INCFD.LT.1 . C IERR = -3 if the X-array is not strictly increasing. C IERR = -4 if IA or IB is out of range. C (Value has not been computed in any of these cases.) C NOTE: The above errors are checked in the order listed, C and following arguments have **NOT** been validated. C C***REFERENCES (NONE) C***ROUTINES CALLED XERROR C***END PROLOGUE DPCHID C C ---------------------------------------------------------------------- C C Change record: C 82-08-04 Converted to SLATEC library version. C 87-07-07 Corrected XERROR calls for d.p. name(s). C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a single precision version, simply: C a. Change DPCHID to PCHID wherever it occurs, C b. Change the double precision declarations to real, and C c. Change the constants ZERO, HALF, SIX to single precision. C C DECLARE ARGUMENTS. C INTEGER N, INCFD, IA, IB, IERR DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N) LOGICAL SKIP C C DECLARE LOCAL VARIABLES. C INTEGER I, IUP, LOW DOUBLE PRECISION H, HALF, SIX, SUM, VALUE, ZERO C C INITIALIZE. C DATA ZERO /0.D0/, HALF/.5D0/, SIX/6.D0/ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT DPCHID IF (SKIP) GO TO 5 C IF ( N.LT.2 ) GO TO 5001 IF ( INCFD.LT.1 ) GO TO 5002 DO 1 I = 2, N IF ( X(I).LE.X(I-1) ) GO TO 5003 1 CONTINUE C C FUNCTION DEFINITION IS OK, GO ON. C 5 CONTINUE SKIP = .TRUE. IF ((IA.LT.1) .OR. (IA.GT.N)) GO TO 5004 IF ((IB.LT.1) .OR. (IB.GT.N)) GO TO 5004 IERR = 0 C C COMPUTE INTEGRAL VALUE. C IF (IA .EQ. IB) THEN VALUE = ZERO ELSE LOW = MIN0(IA, IB) IUP = MAX0(IA, IB) - 1 SUM = ZERO DO 10 I = LOW, IUP H = X(I+1) - X(I) SUM = SUM + H*( (F(1,I) + F(1,I+1)) + * (D(1,I) - D(1,I+1))*(H/SIX) ) 10 CONTINUE VALUE = HALF * SUM IF (IA .GT. IB) VALUE = -VALUE ENDIF C C NORMAL RETURN. C DPCHID = VALUE RETURN C C ERROR RETURNS. C 5001 CONTINUE C N.LT.2 RETURN. IERR = -1 CALL XERROR ('DPCHID -- NUMBER OF DATA POINTS LESS THAN TWO' * , 0, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('DPCHID -- INCREMENT LESS THAN ONE' * , 0, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('DPCHID -- X-ARRAY NOT STRICTLY INCREASING' * , 0, IERR, 1) RETURN C 5004 CONTINUE C IA OR IB OUT OF RANGE RETURN. IERR = -4 CALL XERROR ('DPCHID -- IA OR IB OUT OF RANGE' * , 0, IERR, 1) RETURN C------------- LAST LINE OF DPCHID FOLLOWS ----------------------------- END DOUBLE PRECISION FUNCTION DCHFIV(X1,X2,F1,F2,D1,D2,A,B,IERR) C***BEGIN PROLOGUE DCHFIV C***REFER TO DPCHIA C***ROUTINES CALLED XERROR C***REVISION DATE 870707 (YYMMDD) C***DESCRIPTION C C DCHFIV: Cubic Hermite Function Integral Evaluator. C C Called by DPCHIA to evaluate the integral of a single cubic (in C Hermite form) over an arbitrary interval (A,B). C C ---------------------------------------------------------------------- C C Calling sequence: C C INTEGER IERR C DOUBLE PRECISION X1, X2, F1, F2, D1, D2, A, B C DOUBLE PRECISION VALUE, DCHFIV C C VALUE = DCHFIV (X1, X2, F1, F2, D1, D2, A, B, IERR) C C Parameters: C C VALUE -- (output) VALUE of the requested integral. C C X1,X2 -- (input) endpoints if interval of definition of cubic. C (Must be distinct. Error return if not.) C C F1,F2 -- (input) function values at the ends of the interval. C C D1,D2 -- (input) derivative values at the ends of the interval. C C A,B -- (input) endpoints of interval of integration. C C IERR -- (output) error flag. C Normal return: C IERR = 0 (no errors). C "Recoverable errors": C IERR = -1 if X1.EQ.X2 . C (VALUE has not been set in this case.) C C***END PROLOGUE DCHFIV C C ---------------------------------------------------------------------- C C Programmed by: Fred N. Fritsch, FTS 532-4275, (415) 422-4275, C Mathematics and Statistics Division, C Lawrence Livermore National Laboratory. C C Change record: C 82-08-05 Converted to SLATEC library version. C 87-07-07 Corrected XERROR calls for d.p. name(s). C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a single precision version, simply: C a. Change DCHFIV to CHFIV wherever it occurs, C b. Change the double precision declarations to real, and C c. Change the constants HALF, TWO, ... to single precision. C C DECLARE ARGUMENTS. C INTEGER IERR DOUBLE PRECISION X1, X2, F1, F2, D1, D2, A, B C C DECLARE LOCAL VARIABLES. C DOUBLE PRECISION DTERM, FOUR, FTERM, H, HALF, PHIA1, PHIA2, * PHIB1, PHIB2, PSIA1, PSIA2, PSIB1, PSIB2, TA1, TA2, TB1, * TB2, THREE, TWO, UA1, UA2, UB1, UB2 C C INITIALIZE. C DATA HALF/.5D0/, TWO/2.D0/, THREE/3.D0/, FOUR/4.D0/, SIX/6.D0/ C C VALIDITY CHECK INPUT. C C***FIRST EXECUTABLE STATEMENT DCHFIV IF (X1 .EQ. X2) GO TO 5001 IERR = 0 C C COMPUTE INTEGRAL. C H = X2 - X1 TA1 = (A - X1) / H TA2 = (X2 - A) / H TB1 = (B - X1) / H TB2 = (X2 - B) / H C UA1 = TA1**3 PHIA1 = UA1 * (TWO - TA1) PSIA1 = UA1 * (THREE*TA1 - FOUR) UA2 = TA2**3 PHIA2 = UA2 * (TWO - TA2) PSIA2 = -UA2 * (THREE*TA2 - FOUR) C UB1 = TB1**3 PHIB1 = UB1 * (TWO - TB1) PSIB1 = UB1 * (THREE*TB1 - FOUR) UB2 = TB2**3 PHIB2 = UB2 * (TWO - TB2) PSIB2 = -UB2 * (THREE*TB2 - FOUR) C FTERM = F1*(PHIA2 - PHIB2) + F2*(PHIB1 - PHIA1) DTERM = ( D1*(PSIA2 - PSIB2) + D2*(PSIB1 - PSIA1) )*(H/SIX) C C RETURN VALUE. C DCHFIV = (HALF*H) * (FTERM + DTERM) RETURN C C ERROR RETURN. C 5001 CONTINUE IERR = -1 CALL XERROR ('DCHFIV -- X1 EQUAL TO X2' * , 0, IERR, 1) RETURN C------------- LAST LINE OF DCHFIV FOLLOWS ----------------------------- END