SUBROUTINE DPCHIC(IC,VC,SWITCH,N,X,F,D,INCFD,WK,NWK,IERR) C***BEGIN PROLOGUE DPCHIC C***DATE WRITTEN 820218 (YYMMDD) C***REVISION DATE 870707 (YYMMDD) C***CATEGORY NO. E1B C***KEYWORDS LIBRARY=SLATEC(PCHIP), C TYPE=DOUBLE PRECISION(PCHIC-S DPCHIC-D), C CUBIC HERMITE INTERPOLATION,MONOTONE INTERPOLATION, C PIECEWISE CUBIC INTERPOLATION, C SHAPE-PRESERVING INTERPOLATION C***AUTHOR FRITSCH, F. N., (LLNL) C MATHEMATICS AND STATISTICS DIVISION C LAWRENCE LIVERMORE NATIONAL LABORATORY C P.O. BOX 808 (L-316) C LIVERMORE, CA 94550 C FTS 532-4275, (415) 422-4275 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***DESCRIPTION C C **** Double Precision version of PCHIC **** C C DPCHIC: 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 DPCHFE or DPCHFD. C C ---------------------------------------------------------------------- C C Calling sequence: C C PARAMETER (INCFD = ...) C INTEGER IC(2), N, NWK, IERR C DOUBLE PRECISION VC(2), SWITCH, X(N), F(INCFD,N), D(INCFD,N), C WK(NWK) C C CALL DPCHIC (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 DPCHIM). 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 DPCHSP. 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*8 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 DPCHIM, 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*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). C C D -- (output) real*8 array of derivative values at the data C points. These values will determine a monotone cubic C Hermite function on each subinterval on which the data C are monotonic, except possibly adjacent to switches in C monotonicity. 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*8 array of working storage. The user may C wish 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 AND R.E.CARLSON, 'MONOTONE PIECEWISE C CUBIC INTERPOLATION,' SIAM J.NUMER.ANAL. 17, 2 (APRIL C 1980), 238-246. C 2. F.N.FRITSCH AND J.BUTLAND, 'A METHOD FOR CONSTRUCTING C LOCAL MONOTONE PIECEWISE CUBIC INTERPOLANTS,' SIAM C J.SCI.STAT.COMPUT.5,2 (JUNE 1984), 300-304. C 3. F.N.FRITSCH, 'PIECEWISE CUBIC INTERPOLATION PACKAGE,' C LLNL PREPRINT UCRL-87285 (JULY 1982). C***ROUTINES CALLED DPCHCE,DPCHCI,DPCHCS,XERROR C***END PROLOGUE DPCHIC