LMDZ
chfev.F
Go to the documentation of this file.
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 .LT. 1) GO TO 5001
103  h = x2 - x1
104  IF (h .EQ. 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 500 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.LT.xmi ) next(1) = next(1) + 1
131  IF ( x.GT.xma ) next(2) = next(2) + 1
132 C (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.)
133  500 CONTINUE
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
xermsg
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: xermsg.F:3
chfev
subroutine chfev(X1, X2, F1, F2, D1, D2, NE, XE, FE, NEXT, IERR)
Definition: chfev.F:3