doxy/html/chfev_8_f_source.html

*DECK CHFEV

      SUBROUTINE chfev (X1, X2, F1, F2, D1, D2, NE, XE, FE, NEXT, IERR)

C***BEGIN PROLOGUE  CHFEV

C***PURPOSE  Evaluate a cubic polynomial given in Hermite form at an

C            array of points.  While designed for use by PCHFE, it may

C            be useful directly as an evaluator for a piecewise cubic

C            Hermite function in applications, such as graphing, where

C            the interval is known in advance.

C***LIBRARY   SLATEC (PCHIP)

C***CATEGORY  E3

C***TYPE      SINGLE PRECISION (CHFEV-S, DCHFEV-D)

C***KEYWORDS  CUBIC HERMITE EVALUATION, CUBIC POLYNOMIAL EVALUATION,

C             PCHIP

C***AUTHOR  Fritsch, F. N., (LLNL)

C             Lawrence Livermore National Laboratory

C             P.O. Box 808  (L-316)

C             Livermore, CA  94550

C             FTS 532-4275, (510) 422-4275

C***DESCRIPTION

C

C          CHFEV:  Cubic Hermite Function EValuator

C

C     Evaluates the cubic polynomial determined by function values

C     F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points

C     XE(J), J=1(1)NE.

C

C ----------------------------------------------------------------------

C

C  Calling sequence:

C

C        INTEGER  NE, NEXT(2), IERR

C        REAL  X1, X2, F1, F2, D1, D2, XE(NE), FE(NE)

C

C        CALL  CHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR)

C

C   Parameters:

C

C     X1,X2 -- (input) endpoints of interval of definition of cubic.

C           (Error return if  X1.EQ.X2 .)

C

C     F1,F2 -- (input) values of function at X1 and X2, respectively.

C

C     D1,D2 -- (input) values of derivative at X1 and X2, respectively.

C

C     NE -- (input) number of evaluation points.  (Error return if

C           NE.LT.1 .)

C

C     XE -- (input) real array of points at which the function is to be

C           evaluated.  If any of the XE are outside the interval

C           [X1,X2], a warning error is returned in NEXT.

C

C     FE -- (output) real array of values of the cubic function defined

C           by  X1,X2, F1,F2, D1,D2  at the points  XE.

C

C     NEXT -- (output) integer array indicating number of extrapolation

C           points:

C            NEXT(1) = number of evaluation points to left of interval.

C            NEXT(2) = number of evaluation points to right of interval.

C

C     IERR -- (output) error flag.

C           Normal return:

C              IERR = 0  (no errors).

C           "Recoverable" errors:

C              IERR = -1  if NE.LT.1 .

C              IERR = -2  if X1.EQ.X2 .

C                (The FE-array has not been changed in either case.)

C

C***REFERENCES  (NONE)

C***ROUTINES CALLED  XERMSG

C***REVISION HISTORY  (YYMMDD)

C   811019  DATE WRITTEN

C   820803  Minor cosmetic changes for release 1.

C   890411  Added SAVE statements (Vers. 3.2).

C   890531  Changed all specific intrinsics to generic.  (WRB)

C   890703  Corrected category record.  (WRB)

C   890703  REVISION DATE from Version 3.2

C   891214  Prologue converted to Version 4.0 format.  (BAB)

C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)

C***END PROLOGUE  CHFEV

C  Programming notes:

C

C     To produce a double precision version, simply:

C        a. Change CHFEV to DCHFEV wherever it occurs,

C        b. Change the real declaration to double precision, and

C        c. Change the constant ZERO to double precision.

C

C  DECLARE ARGUMENTS.

C

      INTEGER  ne, next(2), ierr

      REAL  x1, x2, f1, f2, d1, d2, xe(*), fe(*)

C

C  DECLARE LOCAL VARIABLES.

C

      INTEGER  i

      REAL  c2, c3, del1, del2, delta, h, x, xmi, xma, zero

      SAVE zero

      DATA  zero /0./

C

C  VALIDITY-CHECK ARGUMENTS.

C

C***FIRST EXECUTABLE STATEMENT  CHFEV

      IF (ne .LT. 1)  go to 5001

      h = x2 - x1

      IF (h .EQ. zero)  go to 5002

C

C  INITIALIZE.

C

      ierr = 0

      next(1) = 0

      next(2) = 0

      xmi = min(zero, h)

      xma = max(zero, h)

C

C  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1).

C

      delta = (f2 - f1)/h

      del1 = (d1 - delta)/h

      del2 = (d2 - delta)/h

C                                           (DELTA IS NO LONGER NEEDED.)

      c2 = -(del1+del1 + del2)

      c3 = (del1 + del2)/h

C                               (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.)

C

C  EVALUATION LOOP.

C

      DO 500  i = 1, ne

         x = xe(i) - x1

         fe(i) = f1 + x*(d1 + x*(c2 + x*c3))

C          COUNT EXTRAPOLATION POINTS.

         IF ( x.LT.xmi )  next(1) = next(1) + 1

         IF ( x.GT.xma )  next(2) = next(2) + 1

C        (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.)

  500 CONTINUE

C

C  NORMAL RETURN.

C

      RETURN

C

C  ERROR RETURNS.

C

 5001 CONTINUE

C     NE.LT.1 RETURN.

      ierr = -1

      CALL xermsg('SLATEC', 'CHFEV',

     +   'NUMBER OF EVALUATION POINTS LESS THAN ONE', ierr, 1)

      RETURN

C

 5002 CONTINUE

C     X1.EQ.X2 RETURN.

      ierr = -2

      CALL xermsg('SLATEC', 'CHFEV', 'INTERVAL ENDPOINTS EQUAL', ierr,

     +   1)

      RETURN

C------------- LAST LINE OF CHFEV FOLLOWS ------------------------------

      END