!DECK PCHFE SUBROUTINE PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR) !***BEGIN PROLOGUE PCHFE !***PURPOSE Evaluate a piecewise cubic Hermite function at an array of ! points. May be used by itself for Hermite interpolation, ! or as an evaluator for PCHIM or PCHIC. !***LIBRARY SLATEC (PCHIP) !***CATEGORY E3 !***TYPE SINGLE PRECISION (PCHFE-S, DPCHFE-D) !***KEYWORDS CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP, ! PIECEWISE CUBIC EVALUATION !***AUTHOR Fritsch, F. N., (LLNL) ! Lawrence Livermore National Laboratory ! P.O. Box 808 (L-316) ! Livermore, CA 94550 ! FTS 532-4275, (510) 422-4275 !***DESCRIPTION ! ! PCHFE: Piecewise Cubic Hermite Function Evaluator ! ! Evaluates the cubic Hermite function defined by N, X, F, D at ! the points XE(J), J=1(1)NE. ! ! To provide compatibility with PCHIM and PCHIC, includes an ! increment between successive values of the F- and D-arrays. ! ! ---------------------------------------------------------------------- ! ! Calling sequence: ! ! PARAMETER (INCFD = ...) ! INTEGER N, NE, IERR ! REAL X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE) ! LOGICAL SKIP ! ! CALL PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR) ! ! Parameters: ! ! N -- (input) number of data points. (Error return if N.LT.2 .) ! ! X -- (input) real array of independent variable values. The ! elements of X must be strictly increasing: ! X(I-1) .LT. X(I), I = 2(1)N. ! (Error return if not.) ! ! F -- (input) real array of function values. F(1+(I-1)*INCFD) is ! the value corresponding to X(I). ! ! D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is ! the value corresponding to X(I). ! ! INCFD -- (input) increment between successive values in F and D. ! (Error return if INCFD.LT.1 .) ! ! SKIP -- (input/output) logical variable which should be set to ! .TRUE. if the user wishes to skip checks for validity of ! preceding parameters, or to .FALSE. otherwise. ! This will save time in case these checks have already ! been performed (say, in PCHIM or PCHIC). ! SKIP will be set to .TRUE. on normal return. ! ! NE -- (input) number of evaluation points. (Error return if ! NE.LT.1 .) ! ! XE -- (input) real array of points at which the function is to be ! evaluated. ! ! NOTES: ! 1. The evaluation will be most efficient if the elements ! of XE are increasing relative to X; ! that is, XE(J) .GE. X(I) ! implies XE(K) .GE. X(I), all K.GE.J . ! 2. If any of the XE are outside the interval [X(1),X(N)], ! values are extrapolated from the nearest extreme cubic, ! and a warning error is returned. ! ! FE -- (output) real array of values of the cubic Hermite function ! defined by N, X, F, D at the points XE. ! ! IERR -- (output) error flag. ! Normal return: ! IERR = 0 (no errors). ! Warning error: ! IERR.GT.0 means that extrapolation was performed at ! IERR points. ! "Recoverable" errors: ! IERR = -1 if N.LT.2 . ! IERR = -2 if INCFD.LT.1 . ! IERR = -3 if the X-array is not strictly increasing. ! IERR = -4 if NE.LT.1 . ! (The FE-array has not been changed in any of these cases.) ! NOTE: The above errors are checked in the order listed, ! and following arguments have **NOT** been validated. ! !***REFERENCES (NONE) !***ROUTINES CALLED CHFEV, XERMSG !***REVISION HISTORY (YYMMDD) ! 811020 DATE WRITTEN ! 820803 Minor cosmetic changes for release 1. ! 870707 Minor cosmetic changes to prologue. ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) !***END PROLOGUE PCHFE ! Programming notes: ! ! 1. To produce a double precision version, simply: ! a. Change PCHFE to DPCHFE, and CHFEV to DCHFEV, wherever they ! occur, ! b. Change the real declaration to double precision, ! ! 2. Most of the coding between the call to CHFEV and the end of ! the IR-loop could be eliminated if it were permissible to ! assume that XE is ordered relative to X. ! ! 3. CHFEV does not assume that X1 is less than X2. thus, it would ! be possible to write a version of PCHFE that assumes a strict- ! ly decreasing X-array by simply running the IR-loop backwards ! (and reversing the order of appropriate tests). ! ! 4. The present code has a minor bug, which I have decided is not ! worth the effort that would be required to fix it. ! If XE contains points in [X(N-1),X(N)], followed by points .LT. ! X(N-1), followed by points .GT.X(N), the extrapolation points ! will be counted (at least) twice in the total returned in IERR. ! ! DECLARE ARGUMENTS. ! INTEGER :: N, INCFD, NE, IERR REAL :: X(*), F(INCFD,*), D(INCFD,*), XE(*), FE(*) LOGICAL :: SKIP ! ! DECLARE LOCAL VARIABLES. ! INTEGER :: I, IERC, IR, J, JFIRST, NEXT(2), NJ ! ! VALIDITY-CHECK ARGUMENTS. ! !***FIRST EXECUTABLE STATEMENT PCHFE IF (SKIP) GO TO 5 ! IF ( N.LT.2 ) GO TO 5001 IF ( INCFD.LT.1 ) GO TO 5002 DO I = 2, N IF ( X(I).LE.X(I-1) ) GO TO 5003 END DO ! ! FUNCTION DEFINITION IS OK, GO ON. ! 5 CONTINUE IF ( NE.LT.1 ) GO TO 5004 IERR = 0 SKIP = .TRUE. ! ! LOOP OVER INTERVALS. ( INTERVAL INDEX IS IL = IR-1 . ) ! ( INTERVAL IS X(IL).LE.X.LT.X(IR) . ) JFIRST = 1 IR = 2 10 CONTINUE ! ! SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS. ! IF (JFIRST .GT. NE) GO TO 5000 ! ! LOCATE ALL POINTS IN INTERVAL. ! DO J = JFIRST, NE IF (XE(J) .GE. X(IR)) GO TO 30 END DO J = NE + 1 GO TO 40 ! ! HAVE LOCATED FIRST POINT BEYOND INTERVAL. ! 30 CONTINUE IF (IR .EQ. N) J = NE + 1 ! 40 CONTINUE NJ = J - JFIRST ! ! SKIP EVALUATION IF NO POINTS IN INTERVAL. ! IF (NJ .EQ. 0) GO TO 50 ! ! EVALUATE CUBIC AT XE(I), I = JFIRST (1) J-1 . ! ! ---------------------------------------------------------------- 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) ! ---------------------------------------------------------------- IF (IERC .LT. 0) GO TO 5005 ! IF (NEXT(2) .EQ. 0) GO TO 42 ! IF (NEXT(2) .GT. 0) THEN ! IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE ! RIGHT OF X(IR). ! IF (IR .LT. N) GO TO 41 ! IF (IR .EQ. N) THEN ! THESE ARE ACTUALLY EXTRAPOLATION POINTS. IERR = IERR + NEXT(2) GO TO 42 41 CONTINUE ! ELSE ! WE SHOULD NEVER HAVE GOTTEN HERE. GO TO 5005 ! ENDIF ! ENDIF 42 CONTINUE ! IF (NEXT(1) .EQ. 0) GO TO 49 ! IF (NEXT(1) .GT. 0) THEN ! IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE ! LEFT OF X(IR-1). ! IF (IR .GT. 2) GO TO 43 ! IF (IR .EQ. 2) THEN ! THESE ARE ACTUALLY EXTRAPOLATION POINTS. IERR = IERR + NEXT(1) GO TO 49 43 CONTINUE ! ELSE ! XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST ! EVALUATION INTERVAL. ! ! FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1). DO I = JFIRST, J-1 IF (XE(I) .LT. X(IR-1)) GO TO 45 END DO ! NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR ! IN CHFEV. GO TO 5005 ! 45 CONTINUE ! RESET J. (THIS WILL BE THE NEW JFIRST.) J = I ! ! NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY. DO I = 1, IR-1 IF (XE(J) .LT. X(I)) GO TO 47 END DO ! NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1). ! 47 CONTINUE ! AT THIS POINT, EITHER XE(J) .LT. X(1) ! OR X(I-1) .LE. XE(J) .LT. X(I) . ! RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE ! CYCLING. IR = MAX(1, I-1) ! ENDIF ! ENDIF 49 CONTINUE ! JFIRST = J ! ! END OF IR-LOOP. ! 50 CONTINUE IR = IR + 1 IF (IR .LE. N) GO TO 10 ! ! NORMAL RETURN. ! 5000 CONTINUE RETURN ! ! ERROR RETURNS. ! 5001 CONTINUE ! N.LT.2 RETURN. IERR = -1 CALL XERMSG ('SLATEC', 'PCHFE', & 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1) RETURN ! 5002 CONTINUE ! INCFD.LT.1 RETURN. IERR = -2 CALL XERMSG ('SLATEC', 'PCHFE', 'INCREMENT LESS THAN ONE', IERR, & 1) RETURN ! 5003 CONTINUE ! X-ARRAY NOT STRICTLY INCREASING. IERR = -3 CALL XERMSG ('SLATEC', 'PCHFE', 'X-ARRAY NOT STRICTLY INCREASING' & , IERR, 1) RETURN ! 5004 CONTINUE ! NE.LT.1 RETURN. IERR = -4 CALL XERMSG ('SLATEC', 'PCHFE', & 'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1) RETURN ! 5005 CONTINUE ! ERROR RETURN FROM CHFEV. ! *** THIS CASE SHOULD NEVER OCCUR *** IERR = -5 CALL XERMSG ('SLATEC', 'PCHFE', & 'ERROR RETURN FROM CHFEV -- FATAL', IERR, 2) RETURN !------------- LAST LINE OF PCHFE FOLLOWS ------------------------------ END SUBROUTINE PCHFE