1*DECK CHFEV 2 SUBROUTINE CHFEV (X1, X2, F1, F2, D1, D2, NE, XE, FE, NEXT, IERR) 3C***BEGIN PROLOGUE CHFEV 4C***PURPOSE Evaluate a cubic polynomial given in Hermite form at an 5C array of points. While designed for use by PCHFE, 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 SINGLE 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 CHFEV: 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 REAL X1, X2, F1, F2, D1, D2, XE(NE), FE(NE) 33C 34C CALL CHFEV (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 array of points at which the function is to be 49C 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 array of values of the cubic function defined 53C 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 890411 Added SAVE statements (Vers. 3.2). 74C 890531 Changed all specific intrinsics to generic. (WRB) 75C 890703 Corrected category record. (WRB) 76C 890703 REVISION DATE from Version 3.2 77C 891214 Prologue converted to Version 4.0 format. (BAB) 78C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 79C***END PROLOGUE CHFEV 80C Programming notes: 81C 82C To produce a double precision version, simply: 83C a. Change CHFEV to DCHFEV wherever it occurs, 84C b. Change the real declaration to double precision, and 85C c. Change the constant ZERO to double precision. 86C 87C DECLARE ARGUMENTS. 88C 89 INTEGER NE, NEXT(2), IERR 90 REAL X1, X2, F1, F2, D1, D2, XE(*), FE(*) 91C 92C DECLARE LOCAL VARIABLES. 93C 94 INTEGER I 95 REAL C2, C3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO 96 SAVE ZERO 97 DATA ZERO /0./ 98C 99C VALIDITY-CHECK ARGUMENTS. 100C 101C***FIRST EXECUTABLE STATEMENT CHFEV 102 IF (NE .LT. 1) GO TO 5001 103 H = X2 - X1 104 IF (H .EQ. ZERO) GO TO 5002 105C 106C INITIALIZE. 107C 108 IERR = 0 109 NEXT(1) = 0 110 NEXT(2) = 0 111 XMI = MIN(ZERO, H) 112 XMA = MAX(ZERO, H) 113C 114C COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1). 115C 116 DELTA = (F2 - F1)/H 117 DEL1 = (D1 - DELTA)/H 118 DEL2 = (D2 - DELTA)/H 119C (DELTA IS NO LONGER NEEDED.) 120 C2 = -(DEL1+DEL1 + DEL2) 121 C3 = (DEL1 + DEL2)/H 122C (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.) 123C 124C EVALUATION LOOP. 125C 126 DO 500 I = 1, NE 127 X = XE(I) - X1 128 FE(I) = F1 + X*(D1 + X*(C2 + X*C3)) 129C COUNT EXTRAPOLATION POINTS. 130 IF ( X.LT.XMI ) NEXT(1) = NEXT(1) + 1 131 IF ( X.GT.XMA ) NEXT(2) = NEXT(2) + 1 132C (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.) 133 500 CONTINUE 134C 135C NORMAL RETURN. 136C 137 RETURN 138C 139C ERROR RETURNS. 140C 141 5001 CONTINUE 142C NE.LT.1 RETURN. 143 IERR = -1 144 CALL XERMSG ('SLATEC', 'CHFEV', 145 + 'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1) 146 RETURN 147C 148 5002 CONTINUE 149C X1.EQ.X2 RETURN. 150 IERR = -2 151 CALL XERMSG ('SLATEC', 'CHFEV', 'INTERVAL ENDPOINTS EQUAL', IERR, 152 + 1) 153 RETURN 154C------------- LAST LINE OF CHFEV FOLLOWS ------------------------------ 155 END 156