dpchim.f

Go to the documentation of this file.

*DECK DPCHIM
      SUBROUTINE dpchim (N, X, F, D, INCFD, IERR)
C***BEGIN PROLOGUE  DPCHIM
C***PURPOSE  Set derivatives needed to determine a monotone piecewise
C            cubic Hermite interpolant to given data.  Boundary values
C            are provided which are compatible with monotonicity.  The
C            interpolant will have an extremum at each point where mono-
C            tonicity switches direction.  (See DPCHIC if user control
C            is desired over boundary or switch conditions.)
C***LIBRARY   SLATEC (PCHIP)
C***CATEGORY  E1A
C***TYPE      DOUBLE PRECISION (PCHIM-S, DPCHIM-D)
C***KEYWORDS  CUBIC HERMITE INTERPOLATION, MONOTONE INTERPOLATION,
C             PCHIP, PIECEWISE CUBIC INTERPOLATION
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          DPCHIM:  Piecewise Cubic Hermite Interpolation to
C                  Monotone data.
C
C     Sets derivatives needed to determine a monotone piecewise cubic
C     Hermite interpolant to the data given in X and F.
C
C     Default boundary conditions are provided which are compatible
C     with monotonicity.  (See DPCHIC if user control of boundary con-
C     ditions is desired.)
C
C     If the data are only piecewise monotonic, the interpolant will
C     have an extremum at each point where monotonicity switches direc-
C     tion.  (See DPCHIC if user control is desired in such cases.)
C
C     To facilitate two-dimensional applications, includes an increment
C     between successive values of the F- and D-arrays.
C
C     The resulting piecewise cubic Hermite function may be evaluated
C     by DPCHFE or DPCHFD.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  N, IERR
C        DOUBLE PRECISION  X(N), F(INCFD,N), D(INCFD,N)
C
C        CALL  DPCHIM (N, X, F, D, INCFD, IERR)
C
C   Parameters:
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C           If N=2, simply does linear interpolation.
C
C     X -- (input) real*8 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*8 array of dependent variable values to be
C           interpolated.  F(1+(I-1)*INCFD) is value corresponding to
C           X(I).  DPCHIM is designed for monotonic data, but it will
C           work for any F-array.  It will force extrema at points where
C           monotonicity switches direction.  If some other treatment of
C           switch points is desired, DPCHIC should be used instead.
C                                     -----
C     D -- (output) real*8 array of derivative values at the data
C           points.  If the data are monotonic, these values will
C           determine a monotone cubic Hermite function.
C           The value corresponding to X(I) is stored in
C                D(1+(I-1)*INCFD),  I=1(1)N.
C           No other entries in D are changed.
C
C     INCFD -- (input) increment between successive values in F and D.
C           This argument is provided primarily for 2-D applications.
C           (Error return if  INCFD.LT.1 .)
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           Warning error:
C              IERR.GT.0  means that IERR switches in the direction
C                 of monotonicity were detected.
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             (The D-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  1. F. N. Fritsch and J. Butland, A method for construc-
C                 ting local monotone piecewise cubic interpolants, SIAM
C                 Journal on Scientific and Statistical Computing 5, 2
C                 (June 1984), pp. 300-304.
C               2. F. N. Fritsch and R. E. Carlson, Monotone piecewise
C                 cubic interpolation, SIAM Journal on Numerical Ana-
C                 lysis 17, 2 (April 1980), pp. 238-246.
C***ROUTINES CALLED  DPCHST, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   811103  DATE WRITTEN
C   820201  1. Introduced  DPCHST  to reduce possible over/under-
C             flow problems.
C           2. Rearranged derivative formula for same reason.
C   820602  1. Modified end conditions to be continuous functions
C             of data when monotonicity switches in next interval.
C           2. Modified formulas so end conditions are less prone
C             of over/underflow problems.
C   820803  Minor cosmetic changes for release 1.
C   870707  Corrected XERROR calls for d.p. name(s).
C   870813  Updated Reference 1.
C   890206  Corrected XERROR calls.
C   890411  Added SAVE statements (Vers. 3.2).
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890703  Corrected category record.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   891006  Cosmetic changes to prologue.  (WRB)
C   891006  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   920429  Revised format and order of references.  (WRB,FNF)
C***END PROLOGUE  DPCHIM
C  Programming notes:
C
C     1. The function  DPCHST(ARG1,ARG2)  is assumed to return zero if
C        either argument is zero, +1 if they are of the same sign, and
C        -1 if they are of opposite sign.
C     2. To produce a single precision version, simply:
C        a. Change DPCHIM to PCHIM wherever it occurs,
C        b. Change DPCHST to PCHST wherever it occurs,
C        c. Change all references to the Fortran intrinsics to their
C           single precision equivalents,
C        d. Change the double precision declarations to real, and
C        e. Change the constants ZERO and THREE to single precision.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  N, INCFD, IERR
      DOUBLE PRECISION  X(*), F(incfd,*), D(incfd,*)
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I, NLESS1
      DOUBLE PRECISION  DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE,
     *      h1, h2, hsum, hsumt3, three, w1, w2, zero
      SAVE zero, three
      DOUBLE PRECISION  DPCHST
      DATA  zero /0.d0/, three/3.d0/
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  DPCHIM
      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
      ierr = 0
      nless1 = n - 1
      h1 = x(2) - x(1)
      del1 = (f(1,2) - f(1,1))/h1
      dsave = del1
C
C  SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION.
C
      IF (nless1 .GT. 1)  go to 10
      d(1,1) = del1
      d(1,n) = del1
      go to 5000
C
C  NORMAL CASE  (N .GE. 3).
C
   10 CONTINUE
      h2 = x(3) - x(2)
      del2 = (f(1,3) - f(1,2))/h2
C
C  SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
C     SHAPE-PRESERVING.
C
      hsum = h1 + h2
      w1 = (h1 + hsum)/hsum
      w2 = -h1/hsum
      d(1,1) = w1*del1 + w2*del2
      IF ( dpchst(d(1,1),del1) .LE. zero)  THEN
         d(1,1) = zero
      ELSE IF ( dpchst(del1,del2) .LT. zero)  THEN
C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
         dmax = three*del1
         IF (abs(d(1,1)) .GT. abs(dmax))  d(1,1) = dmax
      ENDIF
C
C  LOOP THROUGH INTERIOR POINTS.
C
      DO 50  i = 2, nless1
         IF (i .EQ. 2)  go to 40
C
         h1 = h2
         h2 = x(i+1) - x(i)
         hsum = h1 + h2
         del1 = del2
         del2 = (f(1,i+1) - f(1,i))/h2
   40    CONTINUE
C
C        SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC.
C
         d(1,i) = zero
         IF ( dpchst(del1,del2) .LT. 0.)  go to 42
         IF ( dpchst(del1,del2) .EQ. 0.)  go to 41
         go to 45
C
C        COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY.
C
   41    CONTINUE
         IF (del2 .EQ. zero)  go to 50
         IF ( dpchst(dsave,del2) .LT. zero)  ierr = ierr + 1
         dsave = del2
         go to 50
C
   42    CONTINUE
         ierr = ierr + 1
         dsave = del2
         go to 50
C
C        USE BRODLIE MODIFICATION OF BUTLAND FORMULA.
C
   45    CONTINUE
         hsumt3 = hsum+hsum+hsum
         w1 = (hsum + h1)/hsumt3
         w2 = (hsum + h2)/hsumt3
         dmax = max( abs(del1), abs(del2) )
         dmin = min( abs(del1), abs(del2) )
         drat1 = del1/dmax
         drat2 = del2/dmax
         d(1,i) = dmin/(w1*drat1 + w2*drat2)
C
   50 CONTINUE
C
C  SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
C     SHAPE-PRESERVING.
C
      w1 = -h2/hsum
      w2 = (h2 + hsum)/hsum
      d(1,n) = w1*del1 + w2*del2
      IF ( dpchst(d(1,n),del2) .LE. zero)  THEN
         d(1,n) = zero
      ELSE IF ( dpchst(del1,del2) .LT. zero)  THEN
C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
         dmax = three*del2
         IF (abs(d(1,n)) .GT. abs(dmax))  d(1,n) = dmax
      ENDIF
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', 'DPCHIM',
     +   'NUMBER OF DATA POINTS LESS THAN TWO', ierr, 1)
      RETURN
C
 5002 CONTINUE
C     INCFD.LT.1 RETURN.
      ierr = -2
      CALL xermsg('SLATEC', 'DPCHIM', 'INCREMENT LESS THAN ONE', ierr,
     +   1)
      RETURN
C
 5003 CONTINUE
C     X-ARRAY NOT STRICTLY INCREASING.
      ierr = -3
      CALL xermsg('SLATEC', 'DPCHIM',
     +   'X-ARRAY NOT STRICTLY INCREASING', ierr, 1)
      RETURN
C------------- LAST LINE OF DPCHIM FOLLOWS -----------------------------
      END