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