SUBROUTINE CHFDV(X1,X2,F1,F2,D1,D2,NE,XE,FE,DE,NEXT,IERR) C***BEGIN PROLOGUE CHFDV C***DATE WRITTEN 811019 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E3,H1 C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=SINGLE 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 PCHFD, it may be useful directly as an evaluator for C a piecewise cubic Hermite function in applications, such as C graphing, where the interval is known in advance. C If only function values are required, use CHFEV instead. C***DESCRIPTION C C CHFDV: 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 CHFEV, instead. C C ---------------------------------------------------------------------- C C Calling sequence: C C INTEGER NE, NEXT(2), IERR C REAL X1, X2, F1, F2, D1, D2, XE(NE), FE(NE), DE(NE) C C CALL CHFDV (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 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 array of values of the cubic function defined C by X1,X2, F1,F2, D1,D2 at the points XE. C C DE -- (output) real 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 CHFDV C C ---------------------------------------------------------------------- C C Change record: C 82-08-03 Minor cosmetic changes for release 1. C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a double precision version, simply: C a. Change CHFDV to DCHFDV wherever it occurs, C b. Change the real declaration to double precision, C c. Change the constant ZERO to double precision, and C d. Change the names of the Fortran functions: AMAX1, AMIN1. C C DECLARE ARGUMENTS. C INTEGER NE, NEXT(2), IERR REAL X1, X2, F1, F2, D1, D2, XE(NE), FE(NE), DE(NE) C C DECLARE LOCAL VARIABLES. C INTEGER I REAL C2, C2T2, C3, C3T3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO DATA ZERO /0./ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT CHFDV 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 = AMIN1(ZERO, H) XMA = AMAX1(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 ('CHFDV -- NUMBER OF EVALUATION POINTS LESS THAN ONE' * , 50, IERR, 1) RETURN C 5002 CONTINUE C X1.EQ.X2 RETURN. IERR = -2 CALL XERROR ('CHFDV -- INTERVAL ENDPOINTS EQUAL' * , 33, IERR, 1) RETURN C------------- LAST LINE OF CHFDV FOLLOWS ------------------------------ END SUBROUTINE CHFEV(X1,X2,F1,F2,D1,D2,NE,XE,FE,NEXT,IERR) C***BEGIN PROLOGUE CHFEV C***DATE WRITTEN 811019 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E3,H1 C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=SINGLE PRECISION(CHFEV-S DCHFEV-D), C CUBIC HERMITE EVALUATION,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 at an C array of points. While designed for use by PCHFE, it may C be useful directly as an evaluator for a piecewise cubic C Hermite function in applications, such as graphing, where C the interval is known in advance. C***DESCRIPTION C C CHFEV: Cubic Hermite Function EValuator C C Evaluates the cubic polynomial determined by function values C F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points C XE(J), J=1(1)NE. C C ---------------------------------------------------------------------- C C Calling sequence: C C INTEGER NE, NEXT(2), IERR C REAL X1, X2, F1, F2, D1, D2, XE(NE), FE(NE) C C CALL CHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, 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 array of points at which the function is to be C 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 array of values of the cubic function defined C by X1,X2, F1,F2, D1,D2 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 (The FE-array has not been changed in either case.) C C***REFERENCES (NONE) C***ROUTINES CALLED XERROR C***END PROLOGUE CHFEV C C ---------------------------------------------------------------------- C C Change record: C 82-08-03 Minor cosmetic changes for release 1. C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a double precision version, simply: C a. Change CHFEV to DCHFEV wherever it occurs, C b. Change the real declaration to double precision, C c. Change the constant ZERO to double precision, and C d. Change the names of the Fortran functions: AMAX1, AMIN1. C C DECLARE ARGUMENTS. C INTEGER NE, NEXT(2), IERR REAL X1, X2, F1, F2, D1, D2, XE(NE), FE(NE) C C DECLARE LOCAL VARIABLES. C INTEGER I REAL C2, C3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO DATA ZERO /0./ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT CHFEV 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 = AMIN1(ZERO, H) XMA = AMAX1(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) C3 = (DEL1 + DEL2)/H C (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.) C C EVALUATION LOOP. C DO 500 I = 1, NE X = XE(I) - X1 FE(I) = F1 + X*(D1 + X*(C2 + X*C3)) 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 ('CHFEV -- NUMBER OF EVALUATION POINTS LESS THAN ONE' * , 50, IERR, 1) RETURN C 5002 CONTINUE C X1.EQ.X2 RETURN. IERR = -2 CALL XERROR ('CHFEV -- INTERVAL ENDPOINTS EQUAL' * , 33, IERR, 1) RETURN C------------- LAST LINE OF CHFEV FOLLOWS ------------------------------ END REAL FUNCTION CHFIV(X1,X2,F1,F2,D1,D2,A,B,IERR) C***BEGIN PROLOGUE CHFIV C***REFER TO PCHIA C***ROUTINES CALLED XERROR C***REVISION DATE 870707 (YYMMDD) C***DESCRIPTION C C CHFIV: Cubic Hermite Function Integral Evaluator. C C Called by PCHIA 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 REAL X1, X2, F1, F2, D1, D2, A, B C REAL VALUE, CHFIV C C VALUE = CHFIV (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 CHFIV 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 double precision version, simply: C a. Change CHFIV to DCHFIV wherever it occurs, C b. Change the real declarations to double precision, and C c. Change the constants HALF, TWO, ... to double precision. C C DECLARE ARGUMENTS. C INTEGER IERR REAL X1, X2, F1, F2, D1, D2, A, B C C DECLARE LOCAL VARIABLES. C REAL 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 /0.5/, TWO /2./, THREE /3./, FOUR /4./, SIX /6./ C C VALIDITY CHECK INPUT. C C***FIRST EXECUTABLE STATEMENT CHFIV 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 CHFIV = (HALF*H) * (FTERM + DTERM) RETURN C C ERROR RETURN. C 5001 CONTINUE IERR = -1 CALL XERROR ('CHFIV -- X1 EQUAL TO X2' * , 23, IERR, 1) RETURN C------------- LAST LINE OF CHFIV FOLLOWS ------------------------------ END INTEGER FUNCTION CHFMC(D1,D2,DELTA) C***BEGIN PROLOGUE CHFMC C***REFER TO PCHMC C***ROUTINES CALLED R1MACH C***REVISION DATE 870707 (YYMMDD) C***DESCRIPTION C C CHFMC: Cubic Hermite Function Monotonicity Checker. C C Called by PCHMC to determine the monotonicity properties of the C cubic with boundary derivative values D1,D2 and chord slope DELTA. C C ---------------------------------------------------------------------- C C Calling sequence: C C REAL D1, D2, DELTA C INTEGER ISMON, CHFMC C C ISMON = CHFMC (D1, D2, DELTA) C C Parameters: C C D1,D2 -- (input) derivative values at the ends of an interval. C C DELTA -- (input) data slope over that interval. C C ISMON -- (output) integer function value, indicating the monoto- C nicity of the cubic segment: C ISMON = -1 if function is strictly decreasing; C ISMON = 0 if function is constant; C ISMON = 1 if function is strictly increasing; C ISMON = 2 if function is non-monotonic; C ISMON = 3 if unable to determine. C C Fortran intrinsics used: SIGN. C Other routines used: R1MACH. C C***END PROLOGUE CHFMC 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 83-12-01 Changed from ISIGN to SIGN to correct bug that C produced wrong sign when -1 .LT. DELTA .LT. 0 . C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a double precision version, simply: C a. Change CHFMC to DCHFMC wherever it occurs, C b. Change the real declarations to double precision, and C c. Change the constants ZERO, ONE, ... to double precision. C C DECLARE ARGUMENTS. C REAL D1, D2, DELTA C C DECLARE LOCAL VARIABLES. C INTEGER ISMON, ITRUE REAL A, B, EPS, FOUR, ONE, PHI, TEN, THREE, TWO, ZERO C C INITIALIZE. C DATA ZERO /0./, ONE /1.0/, TWO /2./, THREE /3./, FOUR /4./, 1 TEN /10./ C C MACHINE-DEPENDENT PARAMETER -- SHOULD BE ABOUT 10*UROUND. C***FIRST EXECUTABLE STATEMENT CHFMC EPS = TEN*R1MACH(4) C C MAKE THE CHECK. C IF (DELTA .EQ. ZERO) THEN C CASE OF CONSTANT DATA. IF ((D1.EQ.ZERO) .AND. (D2.EQ.ZERO)) THEN ISMON = 0 ELSE ISMON = 2 ENDIF ELSE C DATA IS NOT CONSTANT -- PICK UP SIGN. ITRUE = SIGN (ONE, DELTA) A = D1/DELTA B = D2/DELTA IF ((A.LT.ZERO) .OR. (B.LT.ZERO)) THEN ISMON = 2 ELSE IF ((A.LE.THREE-EPS) .AND. (B.LE.THREE-EPS)) THEN C INSIDE SQUARE (0,3)X(0,3) IMPLIES OK. ISMON = ITRUE ELSE IF ((A.GT.FOUR+EPS) .AND. (B.GT.FOUR+EPS)) THEN C OUTSIDE SQUARE (0,4)X(0,4) IMPLIES NONMONOTONIC. ISMON = 2 ELSE C MUST CHECK AGAINST BOUNDARY OF ELLIPSE. A = A - TWO B = B - TWO PHI = ((A*A + B*B) + A*B) - THREE IF (PHI .LT. -EPS) THEN ISMON = ITRUE ELSE IF (PHI .GT. EPS) THEN ISMON = 2 ELSE C TO CLOSE TO BOUNDARY TO TELL, C IN THE PRESENCE OF ROUND-OFF ERRORS. ISMON = 3 ENDIF ENDIF ENDIF C C RETURN VALUE. C CHFMC = ISMON RETURN C------------- LAST LINE OF CHFMC FOLLOWS ------------------------------ END SUBROUTINE PCHCE(IC,VC,N,X,H,SLOPE,D,INCFD,IERR) C***BEGIN PROLOGUE PCHCE C***REFER TO PCHIC C***ROUTINES CALLED PCHDF,PCHST,XERROR C***REVISION DATE 870707 (YYMMDD) C***DESCRIPTION C C PCHCE: PCHIC End Derivative Setter. C C Called by PCHIC to set end derivatives as requested by the user. C It must be called after interior derivative values have been set. C ----- C C To facilitate two-dimensional applications, includes an increment C between successive values of the D-array. C C ---------------------------------------------------------------------- C C Calling sequence: C C PARAMETER (INCFD = ...) C INTEGER IC(2), N, IERR C REAL VC(2), X(N), H(N), SLOPE(N), D(INCFD,N) C C CALL PCHCE (IC, VC, N, X, H, SLOPE, D, INCFD, 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 ( see prologue to PCHIC for details. ) C C VC -- (input) real array of length 2 specifying desired boundary C values. VC(1) need be set only if IC(1) = 2 or 3 . C VC(2) need be set only if IC(2) = 2 or 3 . C C N -- (input) number of data points. (assumes N.GE.2) C C X -- (input) real array of independent variable values. (the C elements of X are assumed to be strictly increasing.) C C H -- (input) real array of interval lengths. C SLOPE -- (input) real array of data slopes. C If the data are (X(I),Y(I)), I=1(1)N, then these inputs are: C H(I) = X(I+1)-X(I), C SLOPE(I) = (Y(I+1)-Y(I))/H(I), I=1(1)N-1. C C D -- (input) real array of derivative values at the data points. C The value corresponding to X(I) must be stored in C D(1+(I-1)*INCFD), I=1(1)N. C (output) the value of D at X(1) and/or X(N) is changed, if C necessary, to produce the requested boundary conditions. C no other entries in D are changed. C C INCFD -- (input) increment between successive values in D. C This argument is provided primarily for 2-D applications. C C IERR -- (output) error flag. C Normal return: C IERR = 0 (no errors). C Warning errors: C IERR = 1 if IBEG.LT.0 and D(1) had to be adjusted for C monotonicity. C IERR = 2 if IEND.LT.0 and D(1+(N-1)*INCFD) had to be C adjusted for monotonicity. C IERR = 3 if both of the above are true. C C ------- C WARNING: This routine does no validity-checking of arguments. C ------- C C Fortran intrinsics used: ABS, IABS. C C***END PROLOGUE PCHCE 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 Minor corrections made to prologue.. C C ---------------------------------------------------------------------- C C Programming notes: C C 1. The function PCHST(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 double precision version, simply: C a. Change PCHCE to DPCHCE wherever it occurs, C b. Change PCHDF to DPCHCE wherever it occurs, C c. Change PCHST to DPCHST wherever it occurs, C d. Change all references to the Fortran intrinsics to their C double presision equivalents, C e. Change the real declarations to double precision, and C f. Change the constants ZERO, HALF, ... to double precision. C 3. One could reduce the number of arguments and amount of local C storage, at the expense of reduced code clarity, by passing in C the array WK (rather than splitting it into H and SLOPE) and C increasing its length enough to incorporate STEMP and XTEMP. C 4. The two monotonicity checks only use the sufficient conditions. C Thus, it is possible (but unlikely) for a boundary condition to C be changed, even though the original interpolant was monotonic. C (At least the result is a continuous function of the data.) C C DECLARE ARGUMENTS. C INTEGER IC(2), N, INCFD, IERR REAL VC(2), X(N), H(N), SLOPE(N), D(INCFD,N) C C DELCARE LOCAL VARIABLES. C INTEGER IBEG, IEND, IERF, INDEX, J, K REAL HALF, STEMP(3), THREE, TWO, XTEMP(4), ZERO REAL PCHDF, PCHST C C INITIALIZE. C DATA ZERO /0./, HALF /0.5/, TWO /2./, THREE /3./ C C***FIRST EXECUTABLE STATEMENT PCHCE IBEG = IC(1) IEND = IC(2) IERR = 0 C C SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. C IF ( IABS(IBEG).GT.N ) IBEG = 0 IF ( IABS(IEND).GT.N ) IEND = 0 C C TREAT BEGINNING BOUNDARY CONDITION. C IF (IBEG .EQ. 0) GO TO 2000 K = IABS(IBEG) IF (K .EQ. 1) THEN C BOUNDARY VALUE PROVIDED. D(1,1) = VC(1) ELSE IF (K .EQ. 2) THEN C BOUNDARY SECOND DERIVATIVE PROVIDED. D(1,1) = HALF*( (THREE*SLOPE(1) - D(1,2)) - HALF*VC(1)*H(1) ) ELSE IF (K .LT. 5) THEN C USE K-POINT DERIVATIVE FORMULA. C PICK UP FIRST K POINTS, IN REVERSE ORDER. DO 10 J = 1, K INDEX = K-J+1 C INDEX RUNS FROM K DOWN TO 1. XTEMP(J) = X(INDEX) IF (J .LT. K) STEMP(J) = SLOPE(INDEX-1) 10 CONTINUE C ----------------------------- D(1,1) = PCHDF (K, XTEMP, STEMP, IERF) C ----------------------------- IF (IERF .NE. 0) GO TO 5001 ELSE C USE 'NOT A KNOT' CONDITION. D(1,1) = ( THREE*(H(1)*SLOPE(2) + H(2)*SLOPE(1)) * - TWO*(H(1)+H(2))*D(1,2) - H(1)*D(1,3) ) / H(2) ENDIF C IF (IBEG .GT. 0) GO TO 2000 C C CHECK D(1,1) FOR COMPATIBILITY WITH MONOTONICITY. C IF (SLOPE(1) .EQ. ZERO) THEN IF (D(1,1) .NE. ZERO) THEN D(1,1) = ZERO IERR = IERR + 1 ENDIF ELSE IF ( PCHST(D(1,1),SLOPE(1)) .LT. ZERO) THEN D(1,1) = ZERO IERR = IERR + 1 ELSE IF ( ABS(D(1,1)) .GT. THREE*ABS(SLOPE(1)) ) THEN D(1,1) = THREE*SLOPE(1) IERR = IERR + 1 ENDIF C C TREAT END BOUNDARY CONDITION. C 2000 CONTINUE IF (IEND .EQ. 0) GO TO 5000 K = IABS(IEND) IF (K .EQ. 1) THEN C BOUNDARY VALUE PROVIDED. D(1,N) = VC(2) ELSE IF (K .EQ. 2) THEN C BOUNDARY SECOND DERIVATIVE PROVIDED. D(1,N) = HALF*( (THREE*SLOPE(N-1) - D(1,N-1)) + * HALF*VC(2)*H(N-1) ) ELSE IF (K .LT. 5) THEN C USE K-POINT DERIVATIVE FORMULA. C PICK UP LAST K POINTS. DO 2010 J = 1, K INDEX = N-K+J C INDEX RUNS FROM N+1-K UP TO N. XTEMP(J) = X(INDEX) IF (J .LT. K) STEMP(J) = SLOPE(INDEX) 2010 CONTINUE C ----------------------------- D(1,N) = PCHDF (K, XTEMP, STEMP, IERF) C ----------------------------- IF (IERF .NE. 0) GO TO 5001 ELSE C USE 'NOT A KNOT' CONDITION. D(1,N) = ( THREE*(H(N-1)*SLOPE(N-2) + H(N-2)*SLOPE(N-1)) * - TWO*(H(N-1)+H(N-2))*D(1,N-1) - H(N-1)*D(1,N-2) ) * / H(N-2) ENDIF C IF (IEND .GT. 0) GO TO 5000 C C CHECK D(1,N) FOR COMPATIBILITY WITH MONOTONICITY. C IF (SLOPE(N-1) .EQ. ZERO) THEN IF (D(1,N) .NE. ZERO) THEN D(1,N) = ZERO IERR = IERR + 2 ENDIF ELSE IF ( PCHST(D(1,N),SLOPE(N-1)) .LT. ZERO) THEN D(1,N) = ZERO IERR = IERR + 2 ELSE IF ( ABS(D(1,N)) .GT. THREE*ABS(SLOPE(N-1)) ) THEN D(1,N) = THREE*SLOPE(N-1) IERR = IERR + 2 ENDIF C C NORMAL RETURN. C 5000 CONTINUE RETURN C C ERROR RETURN. C 5001 CONTINUE C ERROR RETURN FROM PCHDF. C *** THIS CASE SHOULD NEVER OCCUR *** IERR = -1 CALL XERROR ('PCHCE -- ERROR RETURN FROM PCHDF' * , 32, IERR, 1) RETURN C------------- LAST LINE OF PCHCE FOLLOWS ------------------------------ END SUBROUTINE PCHCI(N,H,SLOPE,D,INCFD) C***BEGIN PROLOGUE PCHCI C***REFER TO PCHIC C***ROUTINES CALLED PCHST C***REVISION DATE 870707 (YYMMDD) C***DESCRIPTION C C PCHCI: PCHIC Initial Derivative Setter. C C Called by PCHIC to set derivatives needed to determine a monotone C piecewise cubic Hermite interpolant to the data. C C Default boundary conditions are provided which are compatible C with monotonicity. If the data are only piecewise monotonic, the C interpolant will have an extremum at each point where monotonicity C switches direction. C C To facilitate two-dimensional applications, includes an increment C between successive values of the D-array. C C The resulting piecewise cubic Hermite function should be identical C (within roundoff error) to that produced by PCHIM. C C ---------------------------------------------------------------------- C C Calling sequence: C C PARAMETER (INCFD = ...) C INTEGER N C REAL H(N), SLOPE(N), D(INCFD,N) C C CALL PCHCI (N, H, SLOPE, D, INCFD) C C Parameters: C C N -- (input) number of data points. C If N=2, simply does linear interpolation. C C H -- (input) real array of interval lengths. C SLOPE -- (input) real array of data slopes. C If the data are (X(I),Y(I)), I=1(1)N, then these inputs are: C H(I) = X(I+1)-X(I), C SLOPE(I) = (Y(I+1)-Y(I))/H(I), I=1(1)N-1. C C D -- (output) real array of derivative values at the data points. C If the data are monotonic, these values will determine a C 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 D. C This argument is provided primarily for 2-D applications. C C ------- C WARNING: This routine does no validity-checking of arguments. C ------- C C Fortran intrinsics used: ABS, AMAX1, AMIN1. C C***END PROLOGUE PCHCI 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-06-01 Modified end conditions to be continuous functions of C data when monotonicity switches in next interval. C 82-06-02 1. Modified formulas so end conditions are less prone C to over/underflow problems. C 2. Minor modification to HSUM calculation. C 82-08-05 Converted to SLATEC library version. C C ---------------------------------------------------------------------- C C Programming notes: C C 1. The function PCHST(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 double precision version, simply: C a. Change PCHCI to DPCHCI wherever it occurs, C b. Change PCHST to DPCHST wherever it occurs, C c. Change all references to the Fortran intrinsics to their C double presision equivalents, C d. Change the real declarations to double precision, and C e. Change the constants ZERO and THREE to double precision. C C DECLARE ARGUMENTS. C INTEGER N, INCFD REAL H(N), SLOPE(N), D(INCFD,N) C C DECLARE LOCAL VARIBLES. C INTEGER I, NLESS1 REAL DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, HSUM, HSUMT3, THREE, * W1, W2, ZERO REAL PCHST C C INITIALIZE. C DATA ZERO /0./, THREE /3./ C***FIRST EXECUTABLE STATEMENT PCHCI NLESS1 = N - 1 DEL1 = SLOPE(1) 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 DEL2 = SLOPE(2) C C SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE C SHAPE-PRESERVING. C HSUM = H(1) + H(2) W1 = (H(1) + HSUM)/HSUM W2 = -H(1)/HSUM D(1,1) = W1*DEL1 + W2*DEL2 IF ( PCHST(D(1,1),DEL1) .LE. ZERO) THEN D(1,1) = ZERO ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO) THEN C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. DMAX = THREE*DEL1 IF (ABS(D(1,1)) .GT. ABS(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 HSUM = H(I-1) + H(I) DEL1 = DEL2 DEL2 = SLOPE(I) 40 CONTINUE C C SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. C D(1,I) = ZERO IF ( PCHST(DEL1,DEL2) .LE. ZERO) GO TO 50 C C USE BRODLIE MODIFICATION OF BUTLAND FORMULA. C HSUMT3 = HSUM+HSUM+HSUM W1 = (HSUM + H(I-1))/HSUMT3 W2 = (HSUM + H(I) )/HSUMT3 DMAX = AMAX1( ABS(DEL1), ABS(DEL2) ) DMIN = AMIN1( ABS(DEL1), ABS(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 = -H(N-1)/HSUM W2 = (H(N-1) + HSUM)/HSUM D(1,N) = W1*DEL1 + W2*DEL2 IF ( PCHST(D(1,N),DEL2) .LE. ZERO) THEN D(1,N) = ZERO ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO) THEN C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. DMAX = THREE*DEL2 IF (ABS(D(1,N)) .GT. ABS(DMAX)) D(1,N) = DMAX ENDIF C C NORMAL RETURN. C 5000 CONTINUE RETURN C------------- LAST LINE OF PCHCI FOLLOWS ------------------------------ END SUBROUTINE PCHCS(SWITCH,N,H,SLOPE,D,INCFD,IERR) C***BEGIN PROLOGUE PCHCS C***REFER TO PCHIC C***ROUTINES CALLED PCHST,PCHSW C***DESCRIPTION C C PCHCS: PCHIC Monotonicity Switch Derivative Setter. C C Called by PCHIC to adjust the values of D in the vicinity of a C switch in direction of monotonicity, to produce a more "visually C pleasing" curve than that given by PCHIM . C C ---------------------------------------------------------------------- C C Calling sequence: C C PARAMETER (INCFD = ...) C INTEGER N, IERR C REAL SWITCH, H(N), SLOPE(N), D(INCFD,N) C C CALL PCHCS (SWITCH, N, H, SLOPE, D, INCFD, IERR) C C Parameters: C C SWITCH -- (input) indicates the amount of control desired over C local excursions from data. C C N -- (input) number of data points. (assumes N.GT.2 .) C C H -- (input) real array of interval lengths. C SLOPE -- (input) real array of data slopes. C If the data are (X(I),Y(I)), I=1(1)N, then these inputs are: C H(I) = X(I+1)-X(I), C SLOPE(I) = (Y(I+1)-Y(I))/H(I), I=1(1)N-1. C C D -- (input) real array of derivative values at the data points, C as determined by PCHCI. C (output) derivatives in the vicinity of switches in direction C of monotonicity may be adjusted to produce a more "visually C pleasing" curve. 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 D. C This argument is provided primarily for 2-D applications. C C IERR -- (output) error flag. should be zero. C If negative, trouble in PCHSW. (should never happen.) C C ------- C WARNING: This routine does no validity-checking of arguments. C ------- C C Fortran intrinsics used: ABS, AMAX1, AMIN1. C C***END PROLOGUE PCHCS 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-06-17 Redesigned to (1) fix problem with lack of continuity C approaching a flat-topped peak (2) be cleaner and C easier to verify. C Eliminated subroutines PCHSA and PCHSX in the process. C 82-06-22 1. Limited fact to not exceed one, so computed D is a C convex combination of PCHCI value and PCHSD value. C 2. Changed fudge from 1 to 4 (based on experiments). C 82-06-23 Moved PCHSD to an inline function (eliminating MSWTYP). C 82-08-05 Converted to SLATEC library version. C C ---------------------------------------------------------------------- C C Programming notes: C C 1. The function PCHST(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 double precision version, simply: C a. Change PCHCS to DPCHCS wherever it occurs, C b. Change PCHSD to dpchsd wherever it occurs, C c. Change PCHST to DPCHST wherever it occurs, C d. Change PCHSW to DPCHSW wherever it occurs, C e. Change all references to the Fortran intrinsics to their C double precision equivalents, C f. Change the real declarations to double precision, and C g. Change the constants ZERO and ONE to double precision. C C DECLARE ARGUMENTS. C INTEGER N, INCFD, IERR REAL SWITCH, H(N), SLOPE(N), D(INCFD,N) C C DECLARE LOCAL VARIABLES. C INTEGER I, INDX, K, NLESS1 REAL DEL(3), DEXT, DFLOC, DFMX, FACT, FUDGE, ONE, SLMAX, * WTAVE(2), ZERO REAL PCHST C C DEFINE INLINE FUNCTION FOR WEIGHTED AVERAGE OF SLOPES. C REAL PCHSD, S1, S2, H1, H2 PCHSD(S1,S2,H1,H2) = (H2/(H1+H2))*S1 + (H1/(H1+H2))*S2 C C INITIALIZE. C DATA ZERO /0./, ONE /1./ DATA FUDGE /4./ C***FIRST EXECUTABLE STATEMENT PCHCS IERR = 0 NLESS1 = N - 1 C C LOOP OVER SEGMENTS. C DO 900 I = 2, NLESS1 IF ( PCHST(SLOPE(I-1),SLOPE(I)) ) 100, 300, 900 C -------------------------- C 100 CONTINUE C C....... SLOPE SWITCHES MONOTONICITY AT I-TH POINT ..................... C C DO NOT CHANGE D IF 'UP-DOWN-UP'. IF (I .GT. 2) THEN IF ( PCHST(SLOPE(I-2),SLOPE(I)) .GT. ZERO) GO TO 900 C -------------------------- ENDIF IF (I .LT. NLESS1) THEN IF ( PCHST(SLOPE(I+1),SLOPE(I-1)) .GT. ZERO) GO TO 900 C ---------------------------- ENDIF C C ....... COMPUTE PROVISIONAL VALUE FOR D(1,I). C DEXT = PCHSD (SLOPE(I-1), SLOPE(I), H(I-1), H(I)) C C ....... DETERMINE WHICH INTERVAL CONTAINS THE EXTREMUM. C IF ( PCHST(DEXT, SLOPE(I-1)) ) 200, 900, 250 C ----------------------- C 200 CONTINUE C DEXT AND SLOPE(I-1) HAVE OPPOSITE SIGNS -- C EXTREMUM IS IN (X(I-1),X(I)). K = I-1 C SET UP TO COMPUTE NEW VALUES FOR D(1,I-1) AND D(1,I). WTAVE(2) = DEXT IF (K .GT. 1) * WTAVE(1) = PCHSD (SLOPE(K-1), SLOPE(K), H(K-1), H(K)) GO TO 400 C 250 CONTINUE C DEXT AND SLOPE(I) HAVE OPPOSITE SIGNS -- C EXTREMUM IS IN (X(I),X(I+1)). K = I C SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1). WTAVE(1) = DEXT IF (K .LT. NLESS1) * WTAVE(2) = PCHSD (SLOPE(K), SLOPE(K+1), H(K), H(K+1)) GO TO 400 C 300 CONTINUE C C....... AT LEAST ONE OF SLOPE(I-1) AND SLOPE(I) IS ZERO -- C CHECK FOR FLAT-TOPPED PEAK ....................... C IF (I .EQ. NLESS1) GO TO 900 IF ( PCHST(SLOPE(I-1), SLOPE(I+1)) .GE. ZERO) GO TO 900 C ----------------------------- C C WE HAVE FLAT-TOPPED PEAK ON (X(I),X(I+1)). K = I C SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1). WTAVE(1) = PCHSD (SLOPE(K-1), SLOPE(K), H(K-1), H(K)) WTAVE(2) = PCHSD (SLOPE(K), SLOPE(K+1), H(K), H(K+1)) C 400 CONTINUE C C....... AT THIS POINT WE HAVE DETERMINED THAT THERE WILL BE AN EXTREMUM C ON (X(K),X(K+1)), WHERE K=I OR I-1, AND HAVE SET ARRAY WTAVE-- C WTAVE(1) IS A WEIGHTED AVERAGE OF SLOPE(K-1) AND SLOPE(K), C IF K.GT.1 C WTAVE(2) IS A WEIGHTED AVERAGE OF SLOPE(K) AND SLOPE(K+1), C IF K.LT.N-1 C SLMAX = ABS(SLOPE(K)) IF (K .GT. 1) SLMAX = AMAX1( SLMAX, ABS(SLOPE(K-1)) ) IF (K.LT.NLESS1) SLMAX = AMAX1( SLMAX, ABS(SLOPE(K+1)) ) C IF (K .GT. 1) DEL(1) = SLOPE(K-1) / SLMAX DEL(2) = SLOPE(K) / SLMAX IF (K.LT.NLESS1) DEL(3) = SLOPE(K+1) / SLMAX C IF ((K.GT.1) .AND. (K.LT.NLESS1)) THEN C NORMAL CASE -- EXTREMUM IS NOT IN A BOUNDARY INTERVAL. FACT = FUDGE* ABS(DEL(3)*(DEL(1)-DEL(2))*(WTAVE(2)/SLMAX)) D(1,K) = D(1,K) + AMIN1(FACT,ONE)*(WTAVE(1) - D(1,K)) FACT = FUDGE* ABS(DEL(1)*(DEL(3)-DEL(2))*(WTAVE(1)/SLMAX)) D(1,K+1) = D(1,K+1) + AMIN1(FACT,ONE)*(WTAVE(2) - D(1,K+1)) ELSE C SPECIAL CASE K=1 (WHICH CAN OCCUR ONLY IF I=2) OR C K=NLESS1 (WHICH CAN OCCUR ONLY IF I=NLESS1). FACT = FUDGE* ABS(DEL(2)) D(1,I) = AMIN1(FACT,ONE) * WTAVE(I-K+1) C NOTE THAT I-K+1 = 1 IF K=I (=NLESS1), C I-K+1 = 2 IF K=I-1(=1). ENDIF C C C....... ADJUST IF NECESSARY TO LIMIT EXCURSIONS FROM DATA. C IF (SWITCH .LE. ZERO) GO TO 900 C DFLOC = H(K)*ABS(SLOPE(K)) IF (K .GT. 1) DFLOC = AMAX1( DFLOC, H(K-1)*ABS(SLOPE(K-1)) ) IF (K.LT.NLESS1) DFLOC = AMAX1( DFLOC, H(K+1)*ABS(SLOPE(K+1)) ) DFMX = SWITCH*DFLOC INDX = I-K+1 C INDX = 1 IF K=I, 2 IF K=I-1. C --------------------------------------------------------------- CALL PCHSW (DFMX, INDX, D(1,K), D(1,K+1), H(K), SLOPE(K), IERR) C --------------------------------------------------------------- IF (IERR .NE. 0) RETURN C C....... END OF SEGMENT LOOP. C 900 CONTINUE C RETURN C------------- LAST LINE OF PCHCS FOLLOWS ------------------------------ END REAL FUNCTION PCHDF(K,X,S,IERR) C***BEGIN PROLOGUE PCHDF C***REFER TO PCHCE,PCHSP C***ROUTINES CALLED XERROR C***DESCRIPTION C C PCHDF: PCHIP 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 PCHCE and PCHSP 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 PCHDF 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 ---------------------------------------------------------------------- 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 double precision version, simply: C a. Change PCHDF to DPCHDF wherever it occurs, C b. Change the real declarations to double precision, and C c. Change the constant ZERO to double precision. C***END PROLOGUE PCHDF INTEGER K, IERR REAL X(K), S(K) C C DECLARE LOCAL VARIABLES. C INTEGER I, J REAL VALUE, ZERO DATA ZERO /0./ C C CHECK FOR LEGAL VALUE OF K. C C***FIRST EXECUTABLE STATEMENT PCHDF 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 PCHDF = VALUE RETURN C C ERROR RETURN. C 5001 CONTINUE C K.LT.3 RETURN. IERR = -1 CALL XERROR ('PCHDF -- K LESS THAN THREE' * , 26, IERR, 1) PCHDF = ZERO RETURN C------------- LAST LINE OF PCHDF FOLLOWS ------------------------------ END SUBROUTINE PCHEV(N,X,F,D,NVAL,XVAL,FVAL,DVAL,IERR) C***BEGIN PROLOGUE PCHEV 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 XVAL, C easy to use. C***DESCRIPTION C C PCHEV: Piecewise Cubic Hermite or Spline Derivative Evaluator, C 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 XVAL. 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 capabilities. C C ---------------------------------------------------------------------- C C Calling sequence: CALL PCHEV (N, X, F, D, NVAL, XVAL, FVAL, DVAL, IERR) C C INTEGER N, NVAL, IERR C REAL 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) 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. (Error return if not.) C C F -- (input) real array of function values. F(I) is C the value corresponding to X(I). C C D -- (input) real array of derivative values. D(I) is C 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) real array of points at which the functions are to C 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) real array of values of the cubic Hermite function C defined by N, X, F, D at the points XVAL. C C DVAL -- (output) real array of values of the first derivative of C 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 CHFDV. 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 PCHFD C***END PROLOGUE PCHEV INTEGER N, NVAL, IERR REAL 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 PCHEV C CALL PCHFD(N,X,F,D,INCFD,SKIP,NVAL,XVAL,FVAL,DVAL,IERR) C C 5000 CONTINUE RETURN C C------------- LAST LINE OF PCHEV FOLLOWS ------------------------------ END SUBROUTINE PCHEZ(N,X,F,D,SPLINE,WK,LWK,IERR) C***BEGIN PROLOGUE PCHEZ 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 PCHEZ: 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 PCHEZ. Many other choices C are available in the subroutines PCHIC, PCHIM and PCHSP. C C Use PCHEV to evaluate the resulting function and its derivative. C C ---------------------------------------------------------------------- C C Calling sequence: CALL PCHEZ (N, X, F, D, SPLINE, WK, LWK, IERR) C C INTEGER N, IERR, LWK C REAL 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 are C computed from the derivative of a local quadratic unless this C alters monotonicity. C C WK -- (scratch) real work array, which must be declared by the calling 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 that C IERR switches in the direction of monotonicity were detected. C When SPLINE=.FALSE., PCHEZ guarantees that if the input C data is monotone, the interpolant will be too. This warning C is to alert you to the fact that the input data was not 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 PCHIM,PCHSP C***END PROLOGUE PCHEZ INTEGER N, LWK, IERR REAL X(N), F(N), D(N), WK(LWK) LOGICAL SPLINE C C DECLARE LOCAL VARIABLES. C INTEGER IC(2), INCFD REAL VC(2) DATA IC(1) /0/ DATA IC(2) /0/ DATA INCFD /1/ C C C***FIRST EXECUTABLE STATEMENT PCHEZ C IF ( SPLINE ) THEN CALL PCHSP (IC, VC, N, X, F, D, INCFD, WK, LWK, IERR) ELSE CALL PCHIM (N, X, F, D, INCFD, IERR) ENDIF C C ERROR CONDITIONS ALREADY CHECKED IN PCHSP OR PCHIM RETURN C------------- LAST LINE OF PCHEZ FOLLOWS ------------------------------ END SUBROUTINE PCHFD(N,X,F,D,INCFD,SKIP,NE,XE,FE,DE,IERR) C***BEGIN PROLOGUE PCHFD C***DATE WRITTEN 811020 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E3,H1 C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=SINGLE 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 PCHIM C or PCHIC. If only function values are required, use C PCHFE instead. C***DESCRIPTION C C PCHFD: 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 PCHFE, instead. C C To provide compatibility with PCHIM and PCHIC, 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 REAL X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE), DE(NE) C LOGICAL SKIP C C CALL PCHFD (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 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 function values. F(1+(I-1)*INCFD) is C the value corresponding to X(I). C C D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is C 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 PCHIM or PCHIC). 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 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 array of values of the cubic Hermite function C defined by N, X, F, D at the points XE. C C DE -- (output) real 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 CHFDV. NB: this should never happen. C Notify the author **IMMEDIATELY** if it does. C C***REFERENCES (NONE) C***ROUTINES CALLED CHFDV,XERROR C***END PROLOGUE PCHFD C C ---------------------------------------------------------------------- C C Change record: C 82-08-03 Minor cosmetic changes for release 1. C 87-07-07 Minor cosmetic changes to prologue. C C ---------------------------------------------------------------------- C C Programming notes: C C 1. To produce a double precision version, simply: C a. Change PCHFD to DPCHFD, and CHFDV to DCHFDV, wherever they C occur, C b. Change the real declaration to double precision, C C 2. Most of the coding between the call to CHFDV 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. CHFDV does not assume that X1 is less than X2. thus, it would C be possible to write a version of PCHFD 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 REAL 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 PCHFD 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 CHFDV (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 CHFDV. 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 ('PCHFD -- NUMBER OF DATA POINTS LESS THAN TWO' * , 44, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('PCHFD -- INCREMENT LESS THAN ONE' * , 32, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('PCHFD -- X-ARRAY NOT STRICTLY INCREASING' * , 40, IERR, 1) RETURN C 5004 CONTINUE C NE.LT.1 RETURN. IERR = -4 CALL XERROR ('PCHFD -- NUMBER OF EVALUATION POINTS LESS THAN ONE' * , 50, IERR, 1) RETURN C 5005 CONTINUE C ERROR RETURN FROM CHFDV. C *** THIS CASE SHOULD NEVER OCCUR *** IERR = -5 CALL XERROR ('PCHFD -- ERROR RETURN FROM CHFDV -- FATAL' * , 41, IERR, 2) RETURN C------------- LAST LINE OF PCHFD FOLLOWS ------------------------------ END SUBROUTINE PCHFE(N,X,F,D,INCFD,SKIP,NE,XE,FE,IERR) C***BEGIN PROLOGUE PCHFE C***DATE WRITTEN 811020 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E3 C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=SINGLE PRECISION(PCHFE-S DPCHFE-D), C CUBIC HERMITE EVALUATION,HERMITE INTERPOLATION, C 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 at an array of C points. May be used by itself for Hermite interpolation, C or as an evaluator for PCHIM or PCHIC. C***DESCRIPTION C C PCHFE: Piecewise Cubic Hermite Function Evaluator C C Evaluates the cubic Hermite function defined by N, X, F, D at C the points XE(J), J=1(1)NE. C C To provide compatibility with PCHIM and PCHIC, 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 REAL X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE) C LOGICAL SKIP C C CALL PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR) C C Parameters: C C N -- (input) number of data points. (Error return if N.LT.2 .) 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 function values. F(1+(I-1)*INCFD) is C the value corresponding to X(I). C C D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is C 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 PCHIM or PCHIC). 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 array of points at which the function is to be C evaluated. 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 array of values of the cubic Hermite function C defined by N, X, F, D 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 (The FE-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 (NONE) C***ROUTINES CALLED CHFEV,XERROR C***END PROLOGUE PCHFE C C ---------------------------------------------------------------------- C C Change record: C 82-08-03 Minor cosmetic changes for release 1. C 87-07-07 Minor cosmetic changes to prologue. C C ---------------------------------------------------------------------- C C Programming notes: C C 1. To produce a double precision version, simply: C a. Change PCHFE to DPCHFE, and CHFEV to DCHFEV, wherever they C occur, C b. Change the real declaration to double precision, C C 2. Most of the coding between the call to CHFEV 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. CHFEV does not assume that X1 is less than X2. thus, it would C be possible to write a version of PCHFE 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 REAL X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(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 PCHFE 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 CHFEV (X(IR-1),X(IR), F(1,IR-1),F(1,IR), D(1,IR-1),D(1,IR), * NJ, XE(JFIRST), FE(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 CHFEV. 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 ('PCHFE -- NUMBER OF DATA POINTS LESS THAN TWO' * , 44, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('PCHFE -- INCREMENT LESS THAN ONE' * , 32, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('PCHFE -- X-ARRAY NOT STRICTLY INCREASING' * , 40, IERR, 1) RETURN C 5004 CONTINUE C NE.LT.1 RETURN. IERR = -4 CALL XERROR ('PCHFE -- NUMBER OF EVALUATION POINTS LESS THAN ONE' * , 50, IERR, 1) RETURN C 5005 CONTINUE C ERROR RETURN FROM CHFEV. C *** THIS CASE SHOULD NEVER OCCUR *** IERR = -5 CALL XERROR ('PCHFE -- ERROR RETURN FROM CHFEV -- FATAL' * , 41, IERR, 2) RETURN C------------- LAST LINE OF PCHFE FOLLOWS ------------------------------ END REAL FUNCTION PCHIA(N,X,F,D,INCFD,SKIP,A,B,IERR) C***BEGIN PROLOGUE PCHIA C***DATE WRITTEN 820730 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E3,H2A2 C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=SINGLE 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 PCHIA: 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 PCHIM and PCHIC, 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 REAL X(N), F(INCFD,N), D(INCFD,N), A, B C LOGICAL SKIP C C VALUE = PCHIA (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 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 function values. F(1+(I-1)*INCFD) is C the value corresponding to X(I). C C D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is C 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 PCHIM or PCHIC). 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 CHFIV,PCHID,XERROR C***END PROLOGUE PCHIA C C ---------------------------------------------------------------------- C C C Change record: C 82-08-04 Converted to SLATEC library version. C 87-07-07 Corrected double precision conversion instructions. C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a double precision version, simply: C a. Change PCHIA to DPCHIA wherever it occurs, C b. Change PCHID to DPCHID wherever it occurs, C c. Change CHFIV to DCHFIV wherever it occurs, C d. Change the real declarations to double precision, and C e. Change the constant ZERO to double precision. C C DECLARE ARGUMENTS. C INTEGER N, INCFD, IERR REAL 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 REAL VALUE, XA, XB, ZERO REAL CHFIV, PCHID C C INITIALIZE. C DATA ZERO /0./ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT PCHIA 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 = AMIN1 (A, B) XB = AMAX1 (A, B) IF (XB .LE. X(2)) THEN C INTERVAL IS TO LEFT OF X(2), SO USE FIRST CUBIC. C -------------------------------------------- VALUE = CHFIV (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 = CHFIV(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 = CHFIV (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 = PCHID (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 + CHFIV (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 + CHFIV (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 PCHIA = VALUE RETURN C C ERROR RETURNS. C 5001 CONTINUE C N.LT.2 RETURN. IERR = -1 CALL XERROR ('PCHIA -- NUMBER OF DATA POINTS LESS THAN TWO' * , 44, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('PCHIA -- INCREMENT LESS THAN ONE' * , 32, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('PCHIA -- X-ARRAY NOT STRICTLY INCREASING' * , 40, IERR, 1) RETURN C 5004 CONTINUE C TROUBLE IN CHFIV. (SHOULD NEVER OCCUR.) IERR = -4 CALL XERROR ('PCHIA -- TROUBLE IN CHFIV' * , 25, IERR, 1) RETURN C 5005 CONTINUE C TROUBLE IN PCHID. (SHOULD NEVER OCCUR.) IERR = -5 CALL XERROR ('PCHIA -- TROUBLE IN PCHID' * , 25, IERR, 1) RETURN C------------- LAST LINE OF PCHIA FOLLOWS ------------------------------ END SUBROUTINE PCHIC(IC,VC,SWITCH,N,X,F,D,INCFD,WK,NWK,IERR) C***BEGIN PROLOGUE PCHIC C***DATE WRITTEN 820218 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E1B C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=SINGLE PRECISION(PCHIC-S DPCHIC-D), C CUBIC HERMITE INTERPOLATION,MONOTONE INTERPOLATION, C PIECEWISE CUBIC INTERPOLATION, C SHAPE-PRESERVING 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 piecewise monotone C piecewise cubic Hermite interpolant to given data. C User control is available over boundary conditions and/or C treatment of points where monotonicity switches direction. C***DESCRIPTION C C PCHIC: Piecewise Cubic Hermite Interpolation Coefficients. C C Sets derivatives needed to determine a piecewise monotone piece- C wise cubic interpolant to the data given in X and F satisfying the C boundary conditions specified by IC and VC. C C The treatment of points where monotonicity switches direction is C controlled by argument SWITCH. 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 PCHFE or PCHFD. C C ---------------------------------------------------------------------- C C Calling sequence: C C PARAMETER (INCFD = ...) C INTEGER IC(2), N, NWK, IERR C REAL VC(2), SWITCH, X(N), F(INCFD,N), D(INCFD,N), WK(NWK) C C CALL PCHIC (IC, VC, SWITCH, 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 for the default boundary condition (the same as C used by PCHIM). C If IBEG.NE.0, then its sign indicates whether the boundary C derivative is to be adjusted, if necessary, to be C compatible with monotonicity: C IBEG.GT.0 if no adjustment is to be performed. C IBEG.LT.0 if the derivative is to be adjusted for C monotonicity. C C Allowable values for the magnitude of IBEG are: 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 IBEG = 5 to set D(1) so that the second derivative is con- C tinuous at X(2). (Reverts to the default b.c. if N.LT.4.) C This option is somewhat analogous to the "not a knot" C boundary condition provided by PCHSP. C C NOTES (IBEG): C 1. An error return is taken if ABS(IBEG).GT.5 . C 2. Only in case IBEG.LE.0 is it guaranteed that the C interpolant will be monotonic in the first interval. C If the returned value of D(1) lies between zero and C 3*SLOPE(1), the interpolant will be monotonic. This C is **NOT** checked if IBEG.GT.0 . C 3. If IBEG.LT.0 and D(1) had to be changed to achieve mono- C tonicity, a warning error is returned. 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 (IEND): C 1. An error return is taken if ABS(IEND).GT.5 . C 2. Only in case IEND.LE.0 is it guaranteed that the C interpolant will be monotonic in the last interval. C If the returned value of D(1+(N-1)*INCFD) lies between C zero and 3*SLOPE(N-1), the interpolant will be monotonic. C This is **NOT** checked if IEND.GT.0 . C 3. If IEND.LT.0 and D(1+(N-1)*INCFD) had to be changed to C achieve monotonicity, a warning error is returned. C C VC -- (input) real 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 SWITCH -- (input) indicates desired treatment of points where C direction of monotonicity switches: C Set SWITCH to zero if interpolant is required to be mono- C tonic in each interval, regardless of monotonicity of data. C NOTES: C 1. This will cause D to be set to zero at all switch C points, thus forcing extrema there. C 2. The result of using this option with the default boun- C dary conditions will be identical to using PCHIM, but C will generally cost more compute time. C This option is provided only to facilitate comparison C of different switch and/or boundary conditions. C Set SWITCH nonzero to use a formula based on the 3-point C difference formula in the vicinity of switch points. C If SWITCH is positive, the interpolant on each interval C containing an extremum is controlled to not deviate from C the data by more than SWITCH*DFLOC, where DFLOC is the C maximum of the change of F on this interval and its two C immediate neighbors. C If SWITCH is negative, no such control is to be imposed. C C N -- (input) number of data points. (Error return if N.LT.2 .) 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(1+(I-1)*INCFD) is value corresponding to X(I). C C D -- (output) real array of derivative values at the data points. C These values will determine a monotone cubic Hermite func- C tion on each subinterval on which the data are monotonic, C except possibly adjacent to switches in monotonicity. 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 array of working storage. The user may wish C to know that the returned values are: C WK(I) = H(I) = X(I+1) - X(I) ; C WK(N-1+I) = SLOPE(I) = (F(1,I+1) - F(1,I)) / H(I) C for I = 1(1)N-1. C C NWK -- (input) length of work array. C (Error return if NWK.LT.2*(N-1) .) C C IERR -- (output) error flag. C Normal return: C IERR = 0 (no errors). C Warning errors: C IERR = 1 if IBEG.LT.0 and D(1) had to be adjusted for C monotonicity. C IERR = 2 if IEND.LT.0 and D(1+(N-1)*INCFD) had to be C adjusted for monotonicity. C IERR = 3 if both of the above are true. 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 ABS(IBEG).GT.5 . C IERR = -5 if ABS(IEND).GT.5 . C IERR = -6 if both of the above are true. C IERR = -7 if NWK.LT.2*(N-1) . 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 3. F.N.FRITSCH, 'PIECEWISE CUBIC INTERPOLATION PACKAGE,' C LLNL PREPRINT UCRL-87285 (JULY 1982). C***ROUTINES CALLED PCHCE,PCHCI,PCHCS,XERROR C***END PROLOGUE PCHIC C C ---------------------------------------------------------------------- C C Change record: C 82-08-04 Converted to SLATEC library version. C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a double precision version, simply: C a. Change PCHIC to DPCHIC wherever it occurs, C b. Change PCHCE to DPCHCE wherever it occurs, C c. Change PCHCI to DPCHCI wherever it occurs, C d. Change PCHCS to DPCHCS wherever it occurs, C e. Change the real declarations to double precision, and C f. Change the constant ZERO to double precision. C C DECLARE ARGUMENTS. C INTEGER IC(2), N, INCFD, NWK, IERR REAL VC(2), SWITCH, X(N), F(INCFD,N), D(INCFD,N), WK(NWK) C C DECLARE LOCAL VARIABLES. C INTEGER I, IBEG, IEND, NLESS1 REAL ZERO DATA ZERO /0./ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT PCHIC 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 IBEG = IC(1) IEND = IC(2) IERR = 0 IF (IABS(IBEG) .GT. 5) IERR = IERR - 1 IF (IABS(IEND) .GT. 5) IERR = IERR - 2 IF (IERR .LT. 0) GO TO 5004 C C FUNCTION DEFINITION IS OK -- GO ON. C NLESS1 = N - 1 IF ( NWK .LT. 2*NLESS1 ) GO TO 5007 C C SET UP H AND SLOPE ARRAYS. C DO 20 I = 1, NLESS1 WK(I) = X(I+1) - X(I) WK(NLESS1+I) = (F(1,I+1) - F(1,I)) / WK(I) 20 CONTINUE C C SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. C IF (NLESS1 .GT. 1) GO TO 1000 D(1,1) = WK(2) D(1,N) = WK(2) GO TO 3000 C C NORMAL CASE (N .GE. 3) . C 1000 CONTINUE C C SET INTERIOR DERIVATIVES AND DEFAULT END CONDITIONS. C C -------------------------------------- CALL PCHCI (N, WK(1), WK(N), D, INCFD) C -------------------------------------- C C SET DERIVATIVES AT POINTS WHERE MONOTONICITY SWITCHES DIRECTION. C IF (SWITCH .EQ. ZERO) GO TO 3000 C ---------------------------------------------------- CALL PCHCS (SWITCH, N, WK(1), WK(N), D, INCFD, IERR) C ---------------------------------------------------- IF (IERR .NE. 0) GO TO 5008 C C SET END CONDITIONS. C 3000 CONTINUE IF ( (IBEG.EQ.0) .AND. (IEND.EQ.0) ) GO TO 5000 C ------------------------------------------------------- CALL PCHCE (IC, VC, N, X, WK(1), WK(N), D, INCFD, IERR) C ------------------------------------------------------- IF (IERR .LT. 0) GO TO 5009 C C NORMAL RETURN. C 5000 CONTINUE RETURN C C ERROR RETURNS. C 5001 CONTINUE C N.LT.2 RETURN. IERR = -1 CALL XERROR ('PCHIC -- NUMBER OF DATA POINTS LESS THAN TWO' * , 44, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('PCHIC -- INCREMENT LESS THAN ONE' * , 32, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('PCHIC -- X-ARRAY NOT STRICTLY INCREASING' * , 40, IERR, 1) RETURN C 5004 CONTINUE C IC OUT OF RANGE RETURN. IERR = IERR - 3 CALL XERROR ('PCHIC -- IC OUT OF RANGE' * , 24, IERR, 1) RETURN C 5007 CONTINUE C NWK .LT. 2*(N-1) RETURN. IERR = -7 CALL XERROR ('PCHIC -- WORK ARRAY TOO SMALL' * , 29, IERR, 1) RETURN C 5008 CONTINUE C ERROR RETURN FROM PCHCS. IERR = -8 CALL XERROR ('PCHIC -- ERROR RETURN FROM PCHCS' * , 32, IERR, 1) RETURN C 5009 CONTINUE C ERROR RETURN FROM PCHCE. C *** THIS CASE SHOULD NEVER OCCUR *** IERR = -9 CALL XERROR ('PCHIC -- ERROR RETURN FROM PCHCE' * , 32, IERR, 1) RETURN C------------- LAST LINE OF PCHIC FOLLOWS ------------------------------ END REAL FUNCTION PCHID(N,X,F,D,INCFD,SKIP,IA,IB,IERR) C***BEGIN PROLOGUE PCHID C***DATE WRITTEN 820723 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E1B,H2A2 C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=SINGLE 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 C data points. C***DESCRIPTION C C PCHID: 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 PCHIM and PCHIC, 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 REAL X(N), F(INCFD,N), D(INCFD,N) C LOGICAL SKIP C C VALUE = PCHID (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 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 function values. F(1+(I-1)*INCFD) is C the value corresponding to X(I). C C D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is C 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 PCHIM or PCHIC). 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 PCHID C C ---------------------------------------------------------------------- C C Change record: C 82-08-04 Converted to SLATEC library version. C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a double precision version, simply: C a. Change PCHID to DPCHID wherever it occurs, C b. Change the real declarations to double precision, and C c. Change the constants ZERO, HALF, SIX to double precision. C C DECLARE ARGUMENTS. C INTEGER N, INCFD, IA, IB, IERR REAL X(N), F(INCFD,N), D(INCFD,N) LOGICAL SKIP C C DECLARE LOCAL VARIABLES. C INTEGER I, IUP, LOW REAL H, HALF, SIX, SUM, VALUE, ZERO C C INITIALIZE. C DATA ZERO /0./, HALF /0.5/, SIX /6./ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT PCHID 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 PCHID = VALUE RETURN C C ERROR RETURNS. C 5001 CONTINUE C N.LT.2 RETURN. IERR = -1 CALL XERROR ('PCHID -- NUMBER OF DATA POINTS LESS THAN TWO' * , 44, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('PCHID -- INCREMENT LESS THAN ONE' * , 32, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('PCHID -- X-ARRAY NOT STRICTLY INCREASING' * , 40, IERR, 1) RETURN C 5004 CONTINUE C IA OR IB OUT OF RANGE RETURN. IERR = -4 CALL XERROR ('PCHID -- IA OR IB OUT OF RANGE' * , 30, IERR, 1) RETURN C------------- LAST LINE OF PCHID FOLLOWS ------------------------------ END SUBROUTINE PCHIM(N,X,F,D,INCFD,IERR) C***BEGIN PROLOGUE PCHIM C***DATE WRITTEN 811103 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E1B C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=SINGLE PRECISION(PCHIM-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 PCHIC if user control is C desired over boundary or switch conditions.) C***DESCRIPTION C C PCHIM: 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 PCHIC 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 PCHIC 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 PCHFE or PCHFD. C C ---------------------------------------------------------------------- C C Calling sequence: C C PARAMETER (INCFD = ...) C INTEGER N, IERR C REAL X(N), F(INCFD,N), D(INCFD,N) C C CALL PCHIM (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 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(1+(I-1)*INCFD) is value corresponding to X(I). C PCHIM is designed for monotonic data, but it will work for C any F-array. It will force extrema at points where mono- C tonicity switches direction. If some other treatment of C switch points is desired, PCHIC should be used instead. C ----- C D -- (output) real array of derivative values at the data points. C If the data are monotonic, these values will determine a C 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 PCHST,XERROR C***END PROLOGUE PCHIM C C ---------------------------------------------------------------------- C C Change record: C 82-02-01 1. Introduced PCHST 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 C ---------------------------------------------------------------------- C C Programming notes: C C 1. The function PCHST(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 double precision version, simply: C a. Change PCHIM to DPCHIM wherever it occurs, C b. Change PCHST to DPCHST wherever it occurs, C c. Change all references to the Fortran intrinsics to their C double precision equivalents, C d. Change the real declarations to double precision, and C e. Change the constants ZERO and THREE to double precision. C C DECLARE ARGUMENTS. C INTEGER N, INCFD, IERR REAL X(N), F(INCFD,N), D(INCFD,N) C C DECLARE LOCAL VARIABLES. C INTEGER I, NLESS1 REAL DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE, * H1, H2, HSUM, HSUMT3, THREE, W1, W2, ZERO REAL PCHST DATA ZERO /0./, THREE /3./ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT PCHIM 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 ( PCHST(D(1,1),DEL1) .LE. ZERO) THEN D(1,1) = ZERO ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO) THEN C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. DMAX = THREE*DEL1 IF (ABS(D(1,1)) .GT. ABS(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 ( PCHST(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 ( PCHST(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 = AMAX1( ABS(DEL1), ABS(DEL2) ) DMIN = AMIN1( ABS(DEL1), ABS(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 ( PCHST(D(1,N),DEL2) .LE. ZERO) THEN D(1,N) = ZERO ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO) THEN C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. DMAX = THREE*DEL2 IF (ABS(D(1,N)) .GT. ABS(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 ('PCHIM -- NUMBER OF DATA POINTS LESS THAN TWO' * , 44, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('PCHIM -- INCREMENT LESS THAN ONE' * , 32, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('PCHIM -- X-ARRAY NOT STRICTLY INCREASING' * , 40, IERR, 1) RETURN C------------- LAST LINE OF PCHIM FOLLOWS ------------------------------ END SUBROUTINE PCHMC(N,X,F,D,INCFD,SKIP,ISMON,IERR) C***BEGIN PROLOGUE PCHMC C***DATE WRITTEN 820518 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E1B C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=SINGLE PRECISION(PCHMC-S DPCHMC-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 Check a cubic Hermite function for monotonicity. C***DESCRIPTION C C PCHMC: Piecewise Cubic Hermite Monotonicity Checker. C C Checks the cubic Hermite function defined by N, X, F, D for C monotonicity. C C To provide compatibility with PCHIM and PCHIC, 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, ISMON(N), IERR C REAL X(N), F(INCFD,N), D(INCFD,N) C LOGICAL SKIP C C CALL PCHMC (N, X, F, D, INCFD, SKIP, ISMON, IERR) C C Parameters: C C N -- (input) number of data points. (Error return if N.LT.2 .) 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 function values. F(1+(I-1)*INCFD) is C the value corresponding to X(I). C C D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is C 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. C SKIP will be set to .TRUE. on normal return. C C ISMON -- (output) integer array indicating on which intervals the C PCH function defined by N, X, F, D is monotonic. C For data interval [X(I),X(I+1)], C ISMON(I) = -1 if function is strictly decreasing; C ISMON(I) = 0 if function is constant; C ISMON(I) = 1 if function is strictly increasing; C ISMON(I) = 2 if function is non-monotonic; C ISMON(I) = 3 if unable to determine. (This means that C the D-values are near the boundary of the C monotonicity region. A small increase pro- C duces non-monotonicity; decrease, strict C monotonicity.) C The above applies to I=1(1)N-1. ISMON(N) indicates whether C the entire function is monotonic on [X(1),X(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 (The ISMON-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 F.N.FRITSCH AND R.E.CARLSON, 'MONOTONE PIECEWISE CUBIC C INTERPOLATION,' SIAM J.NUMER.ANAL. 17, 2 (APRIL 1980), C 238-246. C***ROUTINES CALLED CHFMC,XERROR C***END PROLOGUE PCHMC C C ---------------------------------------------------------------------- C C Change record: C 82-08-04 Converted to SLATEC library version. C 83-12-01 Reversed order of subscripts of F and D, so that the C routine will work properly when INCFD.GT.1 . C 87-07-07 Minor cosmetic changes to prologue. C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a double precision version, simply: C a. Change PCHMC to DPCHMC wherever it occurs, C b. Change CHFMC to DCHFMC wherever it occurs, and C c. Change the real declarations to double precision. C C DECLARE ARGUMENTS. C INTEGER N, INCFD, ISMON(N), IERR REAL X(N), F(INCFD,N), D(INCFD,N) LOGICAL SKIP C C DECLARE LOCAL VARIABLES. C INTEGER I, NSEG REAL DELTA INTEGER CHFMC C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT PCHMC 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 SKIP = .TRUE. C C FUNCTION DEFINITION IS OK -- GO ON. C 5 CONTINUE NSEG = N - 1 DO 90 I = 1, NSEG DELTA = (F(1,I+1)-F(1,I))/(X(I+1)-X(I)) C ------------------------------- ISMON(I) = CHFMC (D(1,I), D(1,I+1), DELTA) C ------------------------------- IF (I .EQ. 1) THEN ISMON(N) = ISMON(1) ELSE C NEED TO FIGURE OUT CUMULATIVE MONOTONICITY FROM FOLLOWING C 'MULTIPLICATION TABLE'-- C C * I S M O N (I) C * -1 0 1 2 3 C I *--------------------* C S -1 I -1 -1 2 2 3 I C M 0 I -1 0 1 2 3 I C O 1 I 2 1 1 2 3 I C N 2 I 2 2 2 2 2 I C (N) 3 I 3 3 3 2 3 I C *--------------------* C C IF EQUAL OR ALREADY DECLARED NONMONOTONIC, NO CHANGE NEEDED. IF ((ISMON(I).NE.ISMON(N)) .AND. (ISMON(N).NE.2)) THEN IF ( MAX0(ISMON(I), ISMON(N)) .GT. 1) THEN C AT LEAST ONE IS EITHER 'NO' OR 'MAYBE'. IF (ISMON(I) .EQ. 2) THEN ISMON(N) = 2 ELSE ISMON(N) = 3 ENDIF ELSE IF (ISMON(I)*ISMON(N) .LT. 0) THEN C BOTH MONOTONIC, BUT IN OPPOSITE SENSES. ISMON(N) = 2 ELSE C AT THIS POINT, ONE IS ZERO, THE OTHER IS +-1. ISMON(N) = ISMON(N) + ISMON(I) ENDIF ENDIF ENDIF 90 CONTINUE C C NORMAL RETURN. C IERR = 0 RETURN C C ERROR RETURNS. C 5001 CONTINUE C N.LT.2 RETURN. IERR = -1 CALL XERROR ('PCHMC -- NUMBER OF DATA POINTS LESS THAN TWO' * , 44, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('PCHMC -- INCREMENT LESS THAN ONE' * , 32, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('PCHMC -- X-ARRAY NOT STRICTLY INCREASING' * , 40, IERR, 1) RETURN C------------- LAST LINE OF PCHMC FOLLOWS ------------------------------ END REAL FUNCTION PCHQA(N,X,F,D,A,B,IERR) C***BEGIN PROLOGUE PCHQA 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 Hermite C or spline function over an arbitrary interval, easy to use. C***DESCRIPTION C C PCHQA: 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 PCHIA 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 = PCHQA (N, X, F, D, A, B, IERR) C C INTEGER N, IERR C REAL 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) 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 function values. F(I) is C the value corresponding to X(I). C C D -- (input) real 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 PCHIA C***END PROLOGUE PCHQA INTEGER N, IERR REAL X(N), F(N), D(N), A, B C C DECLARE LOCAL VARIABLES. C INTEGER INCFD REAL PCHIA LOGICAL SKIP C C INITIALIZE. C DATA INCFD /1/ DATA SKIP /.TRUE./ C C C***FIRST EXECUTABLE STATEMENT PCHQA PCHQA = PCHIA( N, X, F, D, INCFD, SKIP, A, B, IERR ) C C ERROR MESSAGES ARE FROM LOWER LEVEL ROUTINES RETURN C C------------- LAST LINE OF PCHQA FOLLOWS ------------------------------ END SUBROUTINE PCHSP(IC,VC,N,X,F,D,INCFD,WK,NWK,IERR) C***BEGIN PROLOGUE PCHSP C***DATE WRITTEN 820503 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E1B C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=SINGLE 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 PCHSP: 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 PCHFE or PCHFD. 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 REAL VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK) C C CALL PCHSP (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 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 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(1+(I-1)*INCFD) is value corresponding to X(I). C C D -- (output) real array of derivative values at the data points. C These values will determine the cubic spline interpolant C 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 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 PCHDF,XERROR C***END PROLOGUE PCHSP C C ---------------------------------------------------------------------- C C Change record: C 82-08-04 Converted to SLATEC library version. C 87-07-07 Minor cosmetic changes to prologue. C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a double precision version, simply: C a. Change PCHSP to DPCHSP wherever it occurs, C b. Change the real declarations to double precision, and C c. Change the constants ZERO, HALF, ... to double precision. C C DECLARE ARGUMENTS. C INTEGER IC(2), N, INCFD, NWK, IERR REAL 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 REAL G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), ZERO REAL PCHDF C DATA ZERO /0./, HALF /0.5/, ONE /1./, TWO /2./, THREE /3./ C C VALIDITY-CHECK ARGUMENTS. C C***FIRST EXECUTABLE STATEMENT PCHSP 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) = PCHDF (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) = PCHDF (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 ('PCHSP -- NUMBER OF DATA POINTS LESS THAN TWO' * , 44, IERR, 1) RETURN C 5002 CONTINUE C INCFD.LT.1 RETURN. IERR = -2 CALL XERROR ('PCHSP -- INCREMENT LESS THAN ONE' * , 32, IERR, 1) RETURN C 5003 CONTINUE C X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERROR ('PCHSP -- X-ARRAY NOT STRICTLY INCREASING' * , 40, IERR, 1) RETURN C 5004 CONTINUE C IC OUT OF RANGE RETURN. IERR = IERR - 3 CALL XERROR ('PCHSP -- IC OUT OF RANGE' * , 24, IERR, 1) RETURN C 5007 CONTINUE C NWK TOO SMALL RETURN. IERR = -7 CALL XERROR ('PCHSP -- WORK ARRAY TOO SMALL' * , 29, 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 ('PCHSP -- SINGULAR LINEAR SYSTEM' * , 31, IERR, 1) RETURN C 5009 CONTINUE C ERROR RETURN FROM PCHDF. C *** THIS CASE SHOULD NEVER OCCUR *** IERR = -9 CALL XERROR ('PCHSP -- ERROR RETURN FROM PCHDF' * , 32, IERR, 1) RETURN C------------- LAST LINE OF PCHSP FOLLOWS ------------------------------ END REAL FUNCTION PCHST(ARG1,ARG2) C***BEGIN PROLOGUE PCHST C***REFER TO PCHCE,PCHCI,PCHCS,PCHIM C***ROUTINES CALLED (NONE) C***DESCRIPTION C C PCHST: PCHIP 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 double precision version, simply: C a. Change PCHST to DPCHST wherever it occurs, C b. Change all references to the Fortran intrinsics to their C double presision equivalents, C c. Change the real declarations to double precision, and C d. Change the constants ZERO and ONE to double precision. C***END PROLOGUE PCHST REAL ARG1, ARG2 C C DECLARE LOCAL VARIABLES. C REAL ONE, ZERO DATA ZERO /0./, ONE /1./ C C PERFORM THE TEST. C C***FIRST EXECUTABLE STATEMENT PCHST PCHST = SIGN(ONE,ARG1) * SIGN(ONE,ARG2) IF ((ARG1.EQ.ZERO) .OR. (ARG2.EQ.ZERO)) PCHST = ZERO C RETURN C------------- LAST LINE OF PCHST FOLLOWS ------------------------------ END SUBROUTINE PCHSW(DFMAX,IEXTRM,D1,D2,H,SLOPE,IERR) C***BEGIN PROLOGUE PCHSW C***REFER TO PCHCS C***ROUTINES CALLED R1MACH,XERROR C***REVISION DATE 870707 (YYMMDD) C***DESCRIPTION C C PCHSW: PCHCS Switch Excursion Limiter. C C Called by PCHCS to adjust D1 and D2 if necessary to insure that C the extremum on this interval is not further than DFMAX from the C extreme data value. C C ---------------------------------------------------------------------- C C Calling sequence: C C INTEGER IEXTRM, IERR C REAL DFMAX, D1, D2, H, SLOPE C C CALL PCHSW (DFMAX, IEXTRM, D1, D2, H, SLOPE, IERR) C C Parameters: C C DFMAX -- (input) maximum allowed difference between F(IEXTRM) and C the cubic determined by derivative values D1,D2. (assumes C DFMAX.GT.0.) C C IEXTRM -- (input) index of the extreme data value. (assumes C IEXTRM = 1 or 2 . Any value .NE.1 is treated as 2.) C C D1,D2 -- (input) derivative values at the ends of the interval. C (Assumes D1*D2 .LE. 0.) C (output) may be modified if necessary to meet the restriction C imposed by DFMAX. C C H -- (input) interval length. (Assumes H.GT.0.) C C SLOPE -- (input) data SLOPE on the interval. C C IERR -- (output) error flag. should be zero. C If IERR=-1, assumption on D1 and D2 is not satisfied. C If IERR=-2, quadratic equation locating extremum has C negative descriminant (should never occur). C C ------- C WARNING: This routine does no validity-checking of arguments. C ------- C C C Fortran intrinsics used: ABS, SIGN, SQRT. C C***END PROLOGUE PCHSW 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 Replaced DATA statement for SMALL with a use of R1MACH. C C ---------------------------------------------------------------------- C C Programming notes: C C To produce a double precision version, simply: C a. Change PCHSW to DPCHSW wherever it occurs, C b. Change PCHCS to DPCHCS wherever it occurs, C c. Change R1MACH to D1MACH wherever it occurs, C d. Change all references to the Fortran intrinsics to their C double precision equivalents, C e. Change the real declarations to double precision, and C f. Change constants ZERO, ONE, TWO, ... to double precision. C C DECLARE ARGUMENTS. C INTEGER IEXTRM, IERR REAL DFMAX, D1, D2, H, SLOPE C C DECLARE LOCAL VARIABLES. C REAL CP, DMAX, FACT, LAMBDA, NU, ONE, PHI, RADCAL, RHO, SIGMA, * SMALL, THAT, THIRD, THREE, TWO, ZERO DATA ZERO /0./, ONE /1./, TWO /2./, THREE /3./, FACT /100./ C THIRD SHOULD BE SLIGHTLY LESS THAN 1/3. DATA THIRD /0.33333/ C C NOTATION AND GENERAL REMARKS. C C RHO IS THE RATIO OF THE DATA SLOPE TO THE DERIVATIVE BEING TESTED. C LAMBDA IS THE RATIO OF D2 TO D1. C THAT = T-HAT(RHO) IS THE NORMALIZED LOCATION OF THE EXTREMUM. C PHI IS THE NORMALIZED VALUE OF P(X)-F1 AT X = XHAT = X-HAT(RHO), C WHERE THAT = (XHAT - X1)/H . C THAT IS, P(XHAT)-F1 = D*H*PHI, WHERE D=D1 OR D2. C SIMILARLY, P(XHAT)-F2 = D*H*(PHI-RHO) . C C SMALL SHOULD BE A FEW ORDERS OF MAGNITUDE GREATER THAN MACHEPS. C***FIRST EXECUTABLE STATEMENT PCHSW SMALL = FACT*R1MACH(4) C C DO MAIN CALCULATION. C IF (D1 .EQ. ZERO) THEN C C SPECIAL CASE -- D1.EQ.ZERO . C C IF D2 IS ALSO ZERO, THIS ROUTINE SHOULD NOT HAVE BEEN CALLED. IF (D2 .EQ. ZERO) GO TO 5001 C RHO = SLOPE/D2 C EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 . IF (RHO .GE. THIRD) GO TO 5000 THAT = (TWO*(THREE*RHO-ONE)) / (THREE*(TWO*RHO-ONE)) PHI = THAT**2 * ((THREE*RHO-ONE)/THREE) C C CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 . IF (IEXTRM .NE. 1) PHI = PHI - RHO C C TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY. DMAX = DFMAX / (H*ABS(PHI)) IF (ABS(D2) .GT. DMAX) D2 = SIGN (DMAX, D2) ELSE C RHO = SLOPE/D1 LAMBDA = -D2/D1 IF (D2 .EQ. ZERO) THEN C C SPECIAL CASE -- D2.EQ.ZERO . C C EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 . IF (RHO .GE. THIRD) GO TO 5000 CP = TWO - THREE*RHO NU = ONE - TWO*RHO THAT = ONE / (THREE*NU) ELSE IF (LAMBDA .LE. ZERO) GO TO 5001 C C NORMAL CASE -- D1 AND D2 BOTH NONZERO, OPPOSITE SIGNS. C NU = ONE - LAMBDA - TWO*RHO SIGMA = ONE - RHO CP = NU + SIGMA IF (ABS(NU) .GT. SMALL) THEN RADCAL = (NU - (TWO*RHO+ONE))*NU + SIGMA**2 IF (RADCAL .LT. ZERO) GO TO 5002 THAT = (CP - SQRT(RADCAL)) / (THREE*NU) ELSE THAT = ONE/(TWO*SIGMA) ENDIF ENDIF PHI = THAT*((NU*THAT - CP)*THAT + ONE) C C CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 . IF (IEXTRM .NE. 1) PHI = PHI - RHO C C TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY. DMAX = DFMAX / (H*ABS(PHI)) IF (ABS(D1) .GT. DMAX) THEN D1 = SIGN (DMAX, D1) D2 = -LAMBDA*D1 ENDIF ENDIF C C NORMAL RETURN. C 5000 CONTINUE IERR = 0 RETURN C C ERROR RETURNS. C 5001 CONTINUE C D1 AND D2 BOTH ZERO, OR BOTH NONZERO AND SAME SIGN. IERR = -1 CALL XERROR ('PCHSW -- D1 AND/OR D2 INVALID' * , 29, IERR, 1) RETURN C 5002 CONTINUE C NEGATIVE VALUE OF RADICAL (SHOULD NEVER OCCUR). IERR = -2 CALL XERROR ('PCHSW -- NEGATIVE RADICAL' * , 25, IERR, 1) RETURN C------------- LAST LINE OF PCHSW FOLLOWS ------------------------------ END