1*DECK DCHFEV 2 SUBROUTINE DCHFEV (X1, X2, F1, F2, D1, D2, NE, XE, FE, NEXT, IERR) 3C***BEGIN PROLOGUE DCHFEV 4C***PURPOSE Evaluate a cubic polynomial given in Hermite form at an 5C array of points. While designed for use by DPCHFE, it may 6C be useful directly as an evaluator for a piecewise cubic 7C Hermite function in applications, such as graphing, where 8C the interval is known in advance. 9C***LIBRARY SLATEC (PCHIP) 10C***CATEGORY E3 11C***TYPE DOUBLE PRECISION (CHFEV-S, DCHFEV-D) 12C***KEYWORDS CUBIC HERMITE EVALUATION, CUBIC POLYNOMIAL EVALUATION, 13C PCHIP 14C***AUTHOR Fritsch, F. N., (LLNL) 15C Lawrence Livermore National Laboratory 16C P.O. Box 808 (L-316) 17C Livermore, CA 94550 18C FTS 532-4275, (510) 422-4275 19C***DESCRIPTION 20C 21C DCHFEV: Cubic Hermite Function EValuator 22C 23C Evaluates the cubic polynomial determined by function values 24C F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points 25C XE(J), J=1(1)NE. 26C 27C ---------------------------------------------------------------------- 28C 29C Calling sequence: 30C 31C INTEGER NE, NEXT(2), IERR 32C DOUBLE PRECISION X1, X2, F1, F2, D1, D2, XE(NE), FE(NE) 33C 34C CALL DCHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR) 35C 36C Parameters: 37C 38C X1,X2 -- (input) endpoints of interval of definition of cubic. 39C (Error return if X1.EQ.X2 .) 40C 41C F1,F2 -- (input) values of function at X1 and X2, respectively. 42C 43C D1,D2 -- (input) values of derivative at X1 and X2, respectively. 44C 45C NE -- (input) number of evaluation points. (Error return if 46C NE.LT.1 .) 47C 48C XE -- (input) real*8 array of points at which the function is to 49C be evaluated. If any of the XE are outside the interval 50C [X1,X2], a warning error is returned in NEXT. 51C 52C FE -- (output) real*8 array of values of the cubic function 53C defined by X1,X2, F1,F2, D1,D2 at the points XE. 54C 55C NEXT -- (output) integer array indicating number of extrapolation 56C points: 57C NEXT(1) = number of evaluation points to left of interval. 58C NEXT(2) = number of evaluation points to right of interval. 59C 60C IERR -- (output) error flag. 61C Normal return: 62C IERR = 0 (no errors). 63C "Recoverable" errors: 64C IERR = -1 if NE.LT.1 . 65C IERR = -2 if X1.EQ.X2 . 66C (The FE-array has not been changed in either case.) 67C 68C***REFERENCES (NONE) 69C***ROUTINES CALLED XERMSG 70C***REVISION HISTORY (YYMMDD) 71C 811019 DATE WRITTEN 72C 820803 Minor cosmetic changes for release 1. 73C 870813 Corrected XERROR calls for d.p. names(s). 74C 890206 Corrected XERROR calls. 75C 890411 Added SAVE statements (Vers. 3.2). 76C 890531 Changed all specific intrinsics to generic. (WRB) 77C 890703 Corrected category record. (WRB) 78C 890831 Modified array declarations. (WRB) 79C 891006 Cosmetic changes to prologue. (WRB) 80C 891006 REVISION DATE from Version 3.2 81C 891214 Prologue converted to Version 4.0 format. (BAB) 82C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 83C***END PROLOGUE DCHFEV 84C Programming notes: 85C 86C To produce a single precision version, simply: 87C a. Change DCHFEV to CHFEV wherever it occurs, 88C b. Change the double precision declaration to real, and 89C c. Change the constant ZERO to single precision. 90C 91C DECLARE ARGUMENTS. 92C 93 INTEGER NE, NEXT(2), IERR 94 DOUBLE PRECISION X1, X2, F1, F2, D1, D2, XE(*), FE(*) 95C 96C DECLARE LOCAL VARIABLES. 97C 98 INTEGER I 99 DOUBLE PRECISION C2, C3, DEL1, DEL2, DELTA, H, X, XMI, XMA, 100 * ZERO 101 SAVE ZERO 102 DATA ZERO /0.D0/ 103C 104C VALIDITY-CHECK ARGUMENTS. 105C 106C***FIRST EXECUTABLE STATEMENT DCHFEV 107 IF (NE .LT. 1) GO TO 5001 108 H = X2 - X1 109 IF (H .EQ. ZERO) GO TO 5002 110C 111C INITIALIZE. 112C 113 IERR = 0 114 NEXT(1) = 0 115 NEXT(2) = 0 116 XMI = MIN(ZERO, H) 117 XMA = MAX(ZERO, H) 118C 119C COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1). 120C 121 DELTA = (F2 - F1)/H 122 DEL1 = (D1 - DELTA)/H 123 DEL2 = (D2 - DELTA)/H 124C (DELTA IS NO LONGER NEEDED.) 125 C2 = -(DEL1+DEL1 + DEL2) 126 C3 = (DEL1 + DEL2)/H 127C (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.) 128C 129C EVALUATION LOOP. 130C 131 DO 500 I = 1, NE 132 X = XE(I) - X1 133 FE(I) = F1 + X*(D1 + X*(C2 + X*C3)) 134C COUNT EXTRAPOLATION POINTS. 135 IF ( X.LT.XMI ) NEXT(1) = NEXT(1) + 1 136 IF ( X.GT.XMA ) NEXT(2) = NEXT(2) + 1 137C (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.) 138 500 CONTINUE 139C 140C NORMAL RETURN. 141C 142 RETURN 143C 144C ERROR RETURNS. 145C 146 5001 CONTINUE 147C NE.LT.1 RETURN. 148 IERR = -1 149 CALL XERMSG ('SLATEC', 'DCHFEV', 150 + 'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1) 151 RETURN 152C 153 5002 CONTINUE 154C X1.EQ.X2 RETURN. 155 IERR = -2 156 CALL XERMSG ('SLATEC', 'DCHFEV', 'INTERVAL ENDPOINTS EQUAL', 157 + IERR, 1) 158 RETURN 159C------------- LAST LINE OF DCHFEV FOLLOWS ----------------------------- 160 END 161