*DECK PCHSP

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

C

C

C  FUNCTION DEFINITION IS OK -- GO ON.

C

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,.).

C

C  SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.

C

C

C  SET UP FOR BOUNDARY CONDITIONS.

C

C        PICK UP FIRST IBEG POINTS, IN REVERSE ORDER.

C           INDEX RUNS FROM IBEG DOWN TO 1.

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

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

      ENDIF

C

C        PICK UP LAST IEND POINTS.

C           INDEX RUNS FROM N+1-IEND UP TO N.

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

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

      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

C           NO CONDITION AT LEFT END AND N = 2.

         ELSE

C           NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2.

            D(1,1) =((WK(1,2) + TWO*WK(1,1))*WK(2,2)*WK(1,3)

         ENDIF

C        SLOPE PRESCRIBED AT LEFT END.

      ELSE

C        SECOND DERIVATIVE PRESCRIBED AT LEFT END.

      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

            D(1,J) = G*D(1,J-1)

      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.

C

C           NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2.

C           EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT*

C           NOT-A-KNOT AT LEFT END POINT).

         ELSE

C           NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR  ALSO NOT-A-

C           KNOT AT LEFT END POINT.

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)

         ENDIF

      ELSE

C        SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT.

      ENDIF

C

C  COMPLETE FORWARD PASS OF GAUSS ELIMINATION.

C

C

C  CARRY OUT BACK SUBSTITUTION

C

   30 CONTINUE

C --------------------(  END  CODING FROM CUBSPL )--------------------

C

C  NORMAL RETURN.

C

C

C  ERROR RETURNS.

C

 5001 CONTINUE

C     N.LT.2 RETURN.

      CALL XERMSG ('SLATEC', 'PCHSP',

C

 5002 CONTINUE

C     INCFD.LT.1 RETURN.

      CALL XERMSG ('SLATEC', 'PCHSP', 'INCREMENT LESS THAN ONE', IERR,

C

 5003 CONTINUE

C     X-ARRAY NOT STRICTLY INCREASING.

      CALL XERMSG ('SLATEC', 'PCHSP', 'X-ARRAY NOT STRICTLY INCREASING'

C

 5004 CONTINUE

C     IC OUT OF RANGE RETURN.

C

 5007 CONTINUE

C     NWK TOO SMALL 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). ***

      CALL XERMSG ('SLATEC', 'PCHSP', 'SINGULAR LINEAR SYSTEM', IERR,

C

 5009 CONTINUE

C     ERROR RETURN FROM PCHDF.

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

      CALL XERMSG ('SLATEC', 'PCHSP', 'ERROR RETURN FROM PCHDF', IERR,

C------------- LAST LINE OF PCHSP FOLLOWS ------------------------------

      END