pchsp.F

Go to the documentation of this file.

*DECK PCHSP
      SUBROUTINE pchsp (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
C***BEGIN PROLOGUE  PCHSP
C***PURPOSE  Set derivatives needed to determine the Hermite represen-
C            tation of the cubic spline interpolant to given data, with
C            specified boundary conditions.
C***LIBRARY   SLATEC (PCHIP)
C***CATEGORY  E1A
C***TYPE      SINGLE PRECISION (PCHSP-S, DPCHSP-D)
C***KEYWORDS  CUBIC HERMITE INTERPOLATION, PCHIP,
C             PIECEWISE CUBIC INTERPOLATION, SPLINE 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          PCHSP:   Piecewise Cubic Hermite Spline
C
C     Computes the Hermite representation of the cubic spline inter-
C     polant to the data given in X and F satisfying the boundary
C     conditions specified by IC and VC.
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 PCHFE or PCHFD.
C
C     NOTE:  This is a modified version of C. de Boor's cubic spline
C            routine CUBSPL.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  IC(2), N, NWK, IERR
C        REAL  VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
C
C        CALL  PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
C
C   Parameters:
C
C     IC -- (input) integer array of length 2 specifying desired
C           boundary conditions:
C           IC(1) = IBEG, desired condition at beginning of data.
C           IC(2) = IEND, desired condition at end of data.
C
C           IBEG = 0  to set D(1) so that the third derivative is con-
C              tinuous at X(2).  This is the "not a knot" condition
C              provided by de Boor's cubic spline routine CUBSPL.
C              < This is the default boundary condition. >
C           IBEG = 1  if first derivative at X(1) is given in VC(1).
C           IBEG = 2  if second derivative at X(1) is given in VC(1).
C           IBEG = 3  to use the 3-point difference formula for D(1).
C                     (Reverts to the default b.c. if N.LT.3 .)
C           IBEG = 4  to use the 4-point difference formula for D(1).
C                     (Reverts to the default b.c. if N.LT.4 .)
C          NOTES:
C           1. An error return is taken if IBEG is out of range.
C           2. For the "natural" boundary condition, use IBEG=2 and
C              VC(1)=0.
C
C           IEND may take on the same values as IBEG, but applied to
C           derivative at X(N).  In case IEND = 1 or 2, the value is
C           given in VC(2).
C
C          NOTES:
C           1. An error return is taken if IEND is out of range.
C           2. For the "natural" boundary condition, use IEND=2 and
C              VC(2)=0.
C
C     VC -- (input) real array of length 2 specifying desired boundary
C           values, as indicated above.
C           VC(1) need be set only if IC(1) = 1 or 2 .
C           VC(2) need be set only if IC(2) = 1 or 2 .
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 dependent variable values to be inter-
C           polated.  F(1+(I-1)*INCFD) is value corresponding to X(I).
C
C     D -- (output) real array of derivative values at the data points.
C           These values will determine the cubic spline interpolant
C           with the requested boundary conditions.
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     WK -- (scratch) real array of working storage.
C
C     NWK -- (input) length of work array.
C           (Error return if NWK.LT.2*N .)
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
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 IBEG.LT.0 or IBEG.GT.4 .
C              IERR = -5  if IEND.LT.0 of IEND.GT.4 .
C              IERR = -6  if both of the above are true.
C              IERR = -7  if NWK is too small.
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C             (The D-array has not been changed in any of these cases.)
C              IERR = -8  in case of trouble solving the linear system
C                         for the interior derivative values.
C             (The D-array may have been changed in this case.)
C             (             Do **NOT** use it!                )
C
C***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
C                 Verlag, New York, 1978, pp. 53-59.
C***ROUTINES CALLED  PCHDF, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   820503  DATE WRITTEN
C   820804  Converted to SLATEC library version.
C   870707  Minor cosmetic changes to prologue.
C   890411  Added SAVE statements (Vers. 3.2).
C   890703  Corrected category record.  (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   920429  Revised format and order of references.  (WRB,FNF)
C***END PROLOGUE  PCHSP
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change PCHSP to DPCHSP wherever it occurs,
C        b. Change the real declarations to double precision, and
C        c. Change the constants ZERO, HALF, ... to double precision.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  IC(2), N, INCFD, NWK, IERR
      REAL  VC(2), X(*), F(incfd,*), D(incfd,*), WK(2,*)
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  IBEG, IEND, INDEX, J, NM1
      REAL  G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), ZERO
      SAVE zero, half, one, two, three
      REAL  PCHDF
C
      DATA  zero /0./,  half /0.5/,  one /1./,  two /2./,  three /3./
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  PCHSP
      IF ( n.LT.2 )  GO TO 5001
      IF ( incfd.LT.1 )  GO TO 5002
      DO 1  j = 2, n
         IF ( x(j).LE.x(j-1) )  GO TO 5003
    1 CONTINUE
C
      ibeg = ic(1)
      iend = ic(2)
      ierr = 0
      IF ( (ibeg.LT.0).OR.(ibeg.GT.4) )  ierr = ierr - 1
      IF ( (iend.LT.0).OR.(iend.GT.4) )  ierr = ierr - 2
      IF ( ierr.LT.0 )  GO TO 5004
C
C  FUNCTION DEFINITION IS OK -- GO ON.
C
      IF ( nwk .LT. 2*n )  GO TO 5007
C
C  COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO,
C  COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.).
      DO 5  j=2,n
         wk(1,j) = x(j) - x(j-1)
         wk(2,j) = (f(1,j) - f(1,j-1))/wk(1,j)
    5 CONTINUE
C
C  SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.
C
      IF ( ibeg.GT.n )  ibeg = 0
      IF ( iend.GT.n )  iend = 0
C
C  SET UP FOR BOUNDARY CONDITIONS.
C
      IF ( (ibeg.EQ.1).OR.(ibeg.EQ.2) )  THEN
         d(1,1) = vc(1)
      ELSE IF (ibeg .GT. 2)  THEN
C        PICK UP FIRST IBEG POINTS, IN REVERSE ORDER.
         DO 10  j = 1, ibeg
            index = ibeg-j+1
C           INDEX RUNS FROM IBEG DOWN TO 1.
            xtemp(j) = x(index)
            IF (j .LT. ibeg)  stemp(j) = wk(2,index)
   10    CONTINUE
C                 --------------------------------
         d(1,1) = pchdf(ibeg, xtemp, stemp, ierr)
C                 --------------------------------
         IF (ierr .NE. 0)  GO TO 5009
         ibeg = 1
      ENDIF
C
      IF ( (iend.EQ.1).OR.(iend.EQ.2) )  THEN
         d(1,n) = vc(2)
      ELSE IF (iend .GT. 2)  THEN
C        PICK UP LAST IEND POINTS.
         DO 15  j = 1, iend
            index = n-iend+j
C           INDEX RUNS FROM N+1-IEND UP TO N.
            xtemp(j) = x(index)
            IF (j .LT. iend)  stemp(j) = wk(2,index+1)
   15    CONTINUE
C                 --------------------------------
         d(1,n) = pchdf(iend, xtemp, stemp, ierr)
C                 --------------------------------
         IF (ierr .NE. 0)  GO TO 5009
         iend = 1
      ENDIF
C
C --------------------( BEGIN CODING FROM CUBSPL )--------------------
C
C  **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF
C  F  AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM-
C  INATION, WITH S(J) ENDING UP IN D(1,J), ALL J.
C     WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE.
C
C  CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM
C             WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1)
C
      IF (ibeg .EQ. 0)  THEN
         IF (n .EQ. 2)  THEN
C           NO CONDITION AT LEFT END AND N = 2.
            wk(2,1) = one
            wk(1,1) = one
            d(1,1) = two*wk(2,2)
         ELSE
C           NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2.
            wk(2,1) = wk(1,3)
            wk(1,1) = wk(1,2) + wk(1,3)
            d(1,1) =((wk(1,2) + two*wk(1,1))*wk(2,2)*wk(1,3)
     *                        + wk(1,2)**2*wk(2,3)) / wk(1,1)
         ENDIF
      ELSE IF (ibeg .EQ. 1)  THEN
C        SLOPE PRESCRIBED AT LEFT END.
         wk(2,1) = one
         wk(1,1) = zero
      ELSE
C        SECOND DERIVATIVE PRESCRIBED AT LEFT END.
         wk(2,1) = two
         wk(1,1) = one
         d(1,1) = three*wk(2,2) - half*wk(1,2)*d(1,1)
      ENDIF
C
C  IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND
C  CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH
C  EQUATION READS    WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J).
C
      nm1 = n-1
      IF (nm1 .GT. 1)  THEN
         DO 20 j=2,nm1
            IF (wk(2,j-1) .EQ. zero)  GO TO 5008
            g = -wk(1,j+1)/wk(2,j-1)
            d(1,j) = g*d(1,j-1)
     *                  + three*(wk(1,j)*wk(2,j+1) + wk(1,j+1)*wk(2,j))
            wk(2,j) = g*wk(1,j-1) + two*(wk(1,j) + wk(1,j+1))
   20    CONTINUE
      ENDIF
C
C  CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM
C           (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N)
C
C     IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK-
C     SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT
C     AT THIS POINT.
      IF (iend .EQ. 1)  GO TO 30
C
      IF (iend .EQ. 0)  THEN
         IF (n.EQ.2 .AND. ibeg.EQ.0)  THEN
C           NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2.
            d(1,2) = wk(2,2)
            GO TO 30
         ELSE IF ((n.EQ.2) .OR. (n.EQ.3 .AND. ibeg.EQ.0))  THEN
C           EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT*
C           NOT-A-KNOT AT LEFT END POINT).
            d(1,n) = two*wk(2,n)
            wk(2,n) = one
            IF (wk(2,n-1) .EQ. zero)  GO TO 5008
            g = -one/wk(2,n-1)
         ELSE
C           NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR  ALSO NOT-A-
C           KNOT AT LEFT END POINT.
            g = wk(1,n-1) + wk(1,n)
C           DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES).
            d(1,n) = ((wk(1,n)+two*g)*wk(2,n)*wk(1,n-1)
     *                  + wk(1,n)**2*(f(1,n-1)-f(1,n-2))/wk(1,n-1))/g
            IF (wk(2,n-1) .EQ. zero)  GO TO 5008
            g = -g/wk(2,n-1)
            wk(2,n) = wk(1,n-1)
         ENDIF
      ELSE
C        SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT.
         d(1,n) = three*wk(2,n) + half*wk(1,n)*d(1,n)
         wk(2,n) = two
         IF (wk(2,n-1) .EQ. zero)  GO TO 5008
         g = -one/wk(2,n-1)
      ENDIF
C
C  COMPLETE FORWARD PASS OF GAUSS ELIMINATION.
C
      wk(2,n) = g*wk(1,n-1) + wk(2,n)
      IF (wk(2,n) .EQ. zero)   GO TO 5008
      d(1,n) = (g*d(1,n-1) + d(1,n))/wk(2,n)
C
C  CARRY OUT BACK SUBSTITUTION
C
   30 CONTINUE
      DO 40 j=nm1,1,-1
         IF (wk(2,j) .EQ. zero)  GO TO 5008
         d(1,j) = (d(1,j) - wk(1,j)*d(1,j+1))/wk(2,j)
   40 CONTINUE
C --------------------(  END  CODING FROM CUBSPL )--------------------
C
C  NORMAL RETURN.
C
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     N.LT.2 RETURN.
      ierr = -1
      CALL xermsg ('SLATEC', 'PCHSP',
     +   'NUMBER OF DATA POINTS LESS THAN TWO', ierr, 1)
      RETURN
C
 5002 CONTINUE
C     INCFD.LT.1 RETURN.
      ierr = -2
      CALL xermsg ('SLATEC', 'PCHSP', 'INCREMENT LESS THAN ONE', ierr,
     +   1)
      RETURN
C
 5003 CONTINUE
C     X-ARRAY NOT STRICTLY INCREASING.
      ierr = -3
      CALL xermsg ('SLATEC', 'PCHSP', 'X-ARRAY NOT STRICTLY INCREASING'
     +   , ierr, 1)
      RETURN
C
 5004 CONTINUE
C     IC OUT OF RANGE RETURN.
      ierr = ierr - 3
      CALL xermsg ('SLATEC', 'PCHSP', 'IC OUT OF RANGE', ierr, 1)
      RETURN
C
 5007 CONTINUE
C     NWK TOO SMALL RETURN.
      ierr = -7
      CALL xermsg ('SLATEC', 'PCHSP', 'WORK ARRAY TOO SMALL', ierr, 1)
      RETURN
C
 5008 CONTINUE
C     SINGULAR SYSTEM.
C   *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES   ***
C   *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). ***
      ierr = -8
      CALL xermsg ('SLATEC', 'PCHSP', 'SINGULAR LINEAR SYSTEM', ierr,
     +   1)
      RETURN
C
 5009 CONTINUE
C     ERROR RETURN FROM PCHDF.
C   *** THIS CASE SHOULD NEVER OCCUR ***
      ierr = -9
      CALL xermsg ('SLATEC', 'PCHSP', 'ERROR RETURN FROM PCHDF', ierr,
     +   1)
      RETURN
C------------- LAST LINE OF PCHSP FOLLOWS ------------------------------
      END