doxy/html/pchfe_8_f_source.html

*DECK PCHFE

      SUBROUTINE pchfe (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)

C***BEGIN PROLOGUE  PCHFE

C***PURPOSE  Evaluate a piecewise cubic Hermite function at an array of

C            points.  May be used by itself for Hermite interpolation,

C            or as an evaluator for PCHIM or PCHIC.

C***LIBRARY   SLATEC (PCHIP)

C***CATEGORY  E3

C***TYPE      SINGLE PRECISION (PCHFE-S, DPCHFE-D)

C***KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,

C             PIECEWISE CUBIC EVALUATION

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          PCHFE:  Piecewise Cubic Hermite Function Evaluator

C

C     Evaluates the cubic Hermite function defined by  N, X, F, D  at

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

C

C     To provide compatibility with PCHIM and PCHIC, includes an

C     increment between successive values of the F- and D-arrays.

C

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

C

C  Calling sequence:

C

C        PARAMETER  (INCFD = ...)

C        INTEGER  N, NE, IERR

C        REAL  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE)

C        LOGICAL  SKIP

C

C        CALL  PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)

C

C   Parameters:

C

C     N -- (input) number of data points.  (Error return if N.LT.2 .)

C

C     X -- (input) real array of independent variable values.  The

C           elements of X must be strictly increasing:

C                X(I-1) .LT. X(I),  I = 2(1)N.

C           (Error return if not.)

C

C     F -- (input) real array of function values.  F(1+(I-1)*INCFD) is

C           the value corresponding to X(I).

C

C     D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is

C           the value corresponding to X(I).

C

C     INCFD -- (input) increment between successive values in F and D.

C           (Error return if  INCFD.LT.1 .)

C

C     SKIP -- (input/output) logical variable which should be set to

C           .TRUE. if the user wishes to skip checks for validity of

C           preceding parameters, or to .FALSE. otherwise.

C           This will save time in case these checks have already

C           been performed (say, in PCHIM or PCHIC).

C           SKIP will be set to .TRUE. on normal return.

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.

C

C          NOTES:

C           1. The evaluation will be most efficient if the elements

C              of XE are increasing relative to X;

C              that is,   XE(J) .GE. X(I)

C              implies    XE(K) .GE. X(I),  all K.GE.J .

C           2. If any of the XE are outside the interval [X(1),X(N)],

C              values are extrapolated from the nearest extreme cubic,

C              and a warning error is returned.

C

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

C           defined by  N, X, F, D  at the points  XE.

C

C     IERR -- (output) error flag.

C           Normal return:

C              IERR = 0  (no errors).

C           Warning error:

C              IERR.GT.0  means that extrapolation was performed at

C                 IERR points.

C           "Recoverable" errors:

C              IERR = -1  if N.LT.2 .

C              IERR = -2  if INCFD.LT.1 .

C              IERR = -3  if the X-array is not strictly increasing.

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

C             (The FE-array has not been changed in any of these cases.)

C               NOTE:  The above errors are checked in the order listed,

C                   and following arguments have **NOT** been validated.

C

C***REFERENCES  (NONE)

C***ROUTINES CALLED  CHFEV, XERMSG

C***REVISION HISTORY  (YYMMDD)

C   811020  DATE WRITTEN

C   820803  Minor cosmetic changes for release 1.

C   870707  Minor cosmetic changes to prologue.

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

C   890831  Modified array declarations.  (WRB)

C   890831  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  PCHFE

C  Programming notes:

C

C     1. To produce a double precision version, simply:

C        a. Change PCHFE to DPCHFE, and CHFEV to DCHFEV, wherever they

C           occur,

C        b. Change the real declaration to double precision,

C

C     2. Most of the coding between the call to CHFEV and the end of

C        the IR-loop could be eliminated if it were permissible to

C        assume that XE is ordered relative to X.

C

C     3. CHFEV does not assume that X1 is less than X2.  thus, it would

C        be possible to write a version of PCHFE that assumes a strict-

C        ly decreasing X-array by simply running the IR-loop backwards

C        (and reversing the order of appropriate tests).

C

C     4. The present code has a minor bug, which I have decided is not

C        worth the effort that would be required to fix it.

C        If XE contains points in [X(N-1),X(N)], followed by points .LT.

C        X(N-1), followed by points .GT.X(N), the extrapolation points

C        will be counted (at least) twice in the total returned in IERR.

C

C  DECLARE ARGUMENTS.

C

      INTEGER  n, incfd, ne, ierr

      REAL  x(*), f(incfd,*), d(incfd,*), xe(*), fe(*)

      LOGICAL  skip

C

C  DECLARE LOCAL VARIABLES.

C

      INTEGER  i, ierc, ir, j, jfirst, next(2), nj

C

C  VALIDITY-CHECK ARGUMENTS.

C

C***FIRST EXECUTABLE STATEMENT  PCHFE

      IF (skip)  go to 5

C

      IF ( n.LT.2 )  go to 5001

      IF ( incfd.LT.1 )  go to 5002

      DO 1  i = 2, n

         IF ( x(i).LE.x(i-1) )  go to 5003

    1 CONTINUE

C

C  FUNCTION DEFINITION IS OK, GO ON.

C

    5 CONTINUE

      IF ( ne.LT.1 )  go to 5004

      ierr = 0

      skip = .true.

C

C  LOOP OVER INTERVALS.        (   INTERVAL INDEX IS  IL = IR-1  . )

C                              ( INTERVAL IS X(IL).LE.X.LT.X(IR) . )

      jfirst = 1

      ir = 2

   10 CONTINUE

C

C     SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS.

C

         IF (jfirst .GT. ne)  go to 5000

C

C     LOCATE ALL POINTS IN INTERVAL.

C

         DO 20  j = jfirst, ne

            IF (xe(j) .GE. x(ir))  go to 30

   20    CONTINUE

         j = ne + 1

         go to 40

C

C     HAVE LOCATED FIRST POINT BEYOND INTERVAL.

C

   30    CONTINUE

         IF (ir .EQ. n)  j = ne + 1

C

   40    CONTINUE

         nj = j - jfirst

C

C     SKIP EVALUATION IF NO POINTS IN INTERVAL.

C

         IF (nj .EQ. 0)  go to 50

C

C     EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 .

C

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

        CALL chfev(x(ir-1),x(ir), f(1,ir-1),f(1,ir), d(1,ir-1),d(1,ir),

     *              nj, xe(jfirst), fe(jfirst), next, ierc)

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

         IF (ierc .LT. 0)  go to 5005

C

         IF (next(2) .EQ. 0)  go to 42

C        IF (NEXT(2) .GT. 0)  THEN

C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE

C           RIGHT OF X(IR).

C

            IF (ir .LT. n)  go to 41

C           IF (IR .EQ. N)  THEN

C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.

               ierr = ierr + next(2)

               go to 42

   41       CONTINUE

C           ELSE

C              WE SHOULD NEVER HAVE GOTTEN HERE.

               go to 5005

C           ENDIF

C        ENDIF

   42    CONTINUE

C

         IF (next(1) .EQ. 0)  go to 49

C        IF (NEXT(1) .GT. 0)  THEN

C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE

C           LEFT OF X(IR-1).

C

            IF (ir .GT. 2)  go to 43

C           IF (IR .EQ. 2)  THEN

C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.

               ierr = ierr + next(1)

               go to 49

   43       CONTINUE

C           ELSE

C              XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST

C              EVALUATION INTERVAL.

C

C              FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1).

               DO 44  i = jfirst, j-1

                  IF (xe(i) .LT. x(ir-1))  go to 45

   44          CONTINUE

C              NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR

C                     IN CHFEV.

               go to 5005

C

   45          CONTINUE

C              RESET J.  (THIS WILL BE THE NEW JFIRST.)

               j = i

C

C              NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY.

               DO 46  i = 1, ir-1

                  IF (xe(j) .LT. x(i)) go to 47

   46          CONTINUE

C              NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1).

C

   47          CONTINUE

C              AT THIS POINT, EITHER  XE(J) .LT. X(1)

C                 OR      X(I-1) .LE. XE(J) .LT. X(I) .

C              RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE

C              CYCLING.

               ir = max(1, i-1)

C           ENDIF

C        ENDIF

   49    CONTINUE

C

         jfirst = j

C

C     END OF IR-LOOP.

C

   50 CONTINUE

      ir = ir + 1

      IF (ir .LE. n)  go to 10

C

C  NORMAL RETURN.

C

 5000 CONTINUE

      RETURN

C

C  ERROR RETURNS.

C

 5001 CONTINUE

C     N.LT.2 RETURN.

      ierr = -1

      CALL xermsg('SLATEC', 'PCHFE',

     +   'NUMBER OF DATA POINTS LESS THAN TWO', ierr, 1)

      RETURN

C

 5002 CONTINUE

C     INCFD.LT.1 RETURN.

      ierr = -2

      CALL xermsg('SLATEC', 'PCHFE', 'INCREMENT LESS THAN ONE', ierr,

     +   1)

      RETURN

C

 5003 CONTINUE

C     X-ARRAY NOT STRICTLY INCREASING.

      ierr = -3

      CALL xermsg('SLATEC', 'PCHFE', 'X-ARRAY NOT STRICTLY INCREASING'

     +   , ierr, 1)

      RETURN

C

 5004 CONTINUE

C     NE.LT.1 RETURN.

      ierr = -4

      CALL xermsg('SLATEC', 'PCHFE',

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

      RETURN

C

 5005 CONTINUE

C     ERROR RETURN FROM CHFEV.

C   *** THIS CASE SHOULD NEVER OCCUR ***

      ierr = -5

      CALL xermsg('SLATEC', 'PCHFE',

     +   'ERROR RETURN FROM CHFEV -- FATAL', ierr, 2)

      RETURN

C------------- LAST LINE OF PCHFE FOLLOWS ------------------------------

      END