SLATEC Routines --- PCHIC ---
*DECK PCHIC
SUBROUTINE PCHIC (IC, VC, SWITCH, N, X, F, D, INCFD, WK, NWK,
+ IERR)
C***BEGIN PROLOGUE PCHIC
C***PURPOSE Set derivatives needed to determine a piecewise monotone
C piecewise cubic Hermite interpolant to given data.
C User control is available over boundary conditions and/or
C treatment of points where monotonicity switches direction.
C***LIBRARY SLATEC (PCHIP)
C***CATEGORY E1A
C***TYPE SINGLE PRECISION (PCHIC-S, DPCHIC-D)
C***KEYWORDS CUBIC HERMITE INTERPOLATION, MONOTONE INTERPOLATION,
C PCHIP, PIECEWISE CUBIC INTERPOLATION,
C SHAPE-PRESERVING 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 PCHIC: Piecewise Cubic Hermite Interpolation Coefficients.
C
C Sets derivatives needed to determine a piecewise monotone piece-
C wise cubic interpolant to the data given in X and F satisfying the
C boundary conditions specified by IC and VC.
C
C The treatment of points where monotonicity switches direction is
C controlled by argument SWITCH.
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 ----------------------------------------------------------------------
C
C Calling sequence:
C
C PARAMETER (INCFD = ...)
C INTEGER IC(2), N, NWK, IERR
C REAL VC(2), SWITCH, X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
C
C CALL PCHIC (IC, VC, SWITCH, 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 for the default boundary condition (the same as
C used by PCHIM).
C If IBEG.NE.0, then its sign indicates whether the boundary
C derivative is to be adjusted, if necessary, to be
C compatible with monotonicity:
C IBEG.GT.0 if no adjustment is to be performed.
C IBEG.LT.0 if the derivative is to be adjusted for
C monotonicity.
C
C Allowable values for the magnitude of IBEG are:
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 IBEG = 5 to set D(1) so that the second derivative is con-
C tinuous at X(2). (Reverts to the default b.c. if N.LT.4.)
C This option is somewhat analogous to the "not a knot"
C boundary condition provided by PCHSP.
C
C NOTES (IBEG):
C 1. An error return is taken if ABS(IBEG).GT.5 .
C 2. Only in case IBEG.LE.0 is it guaranteed that the
C interpolant will be monotonic in the first interval.
C If the returned value of D(1) lies between zero and
C 3*SLOPE(1), the interpolant will be monotonic. This
C is **NOT** checked if IBEG.GT.0 .
C 3. If IBEG.LT.0 and D(1) had to be changed to achieve mono-
C tonicity, a warning error is returned.
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 (IEND):
C 1. An error return is taken if ABS(IEND).GT.5 .
C 2. Only in case IEND.LE.0 is it guaranteed that the
C interpolant will be monotonic in the last interval.
C If the returned value of D(1+(N-1)*INCFD) lies between
C zero and 3*SLOPE(N-1), the interpolant will be monotonic.
C This is **NOT** checked if IEND.GT.0 .
C 3. If IEND.LT.0 and D(1+(N-1)*INCFD) had to be changed to
C achieve monotonicity, a warning error is returned.
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 SWITCH -- (input) indicates desired treatment of points where
C direction of monotonicity switches:
C Set SWITCH to zero if interpolant is required to be mono-
C tonic in each interval, regardless of monotonicity of data.
C NOTES:
C 1. This will cause D to be set to zero at all switch
C points, thus forcing extrema there.
C 2. The result of using this option with the default boun-
C dary conditions will be identical to using PCHIM, but
C will generally cost more compute time.
C This option is provided only to facilitate comparison
C of different switch and/or boundary conditions.
C Set SWITCH nonzero to use a formula based on the 3-point
C difference formula in the vicinity of switch points.
C If SWITCH is positive, the interpolant on each interval
C containing an extremum is controlled to not deviate from
C the data by more than SWITCH*DFLOC, where DFLOC is the
C maximum of the change of F on this interval and its two
C immediate neighbors.
C If SWITCH is negative, no such control is to be imposed.
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 a monotone cubic Hermite func-
C tion on each subinterval on which the data are monotonic,
C except possibly adjacent to switches in monotonicity.
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. The user may wish
C to know that the returned values are:
C WK(I) = H(I) = X(I+1) - X(I) ;
C WK(N-1+I) = SLOPE(I) = (F(1,I+1) - F(1,I)) / H(I)
C for I = 1(1)N-1.
C
C NWK -- (input) length of work array.
C (Error return if NWK.LT.2*(N-1) .)
C
C IERR -- (output) error flag.
C Normal return:
C IERR = 0 (no errors).
C Warning errors:
C IERR = 1 if IBEG.LT.0 and D(1) had to be adjusted for
C monotonicity.
C IERR = 2 if IEND.LT.0 and D(1+(N-1)*INCFD) had to be
C adjusted for monotonicity.
C IERR = 3 if both of the above are true.
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 ABS(IBEG).GT.5 .
C IERR = -5 if ABS(IEND).GT.5 .
C IERR = -6 if both of the above are true.
C IERR = -7 if NWK.LT.2*(N-1) .
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, Piecewise Cubic Hermite Interpolation
C Package, Report UCRL-87285, Lawrence Livermore Nation-
C al Laboratory, July 1982. [Poster presented at the
C SIAM 30th Anniversary Meeting, 19-23 July 1982.]
C 2. 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 3. 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 PCHCE, PCHCI, PCHCS, XERMSG
C***REVISION HISTORY (YYMMDD)
C 820218 DATE WRITTEN
C 820804 Converted to SLATEC library version.
C 870813 Updated Reference 2.
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 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 PCHIC