chfev.f

Go to the documentation of this file.

C                          $Author: ulammers $
C                          $Date: 1997/02/14 14:09:43 $
C                          $Id: chfev.f,v 1.1 1997/02/14 14:09:43 ulammers Exp $
C                          $Source: /usr4/users/aparmar/SAXDAS/pipeline/LE_lemat/RCS/chfev.f,v $
*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