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