1*DECK DPCHIM 2 SUBROUTINE DPCHIM (N, X, F, D, INCFD, IERR) 3C***BEGIN PROLOGUE DPCHIM 4C***PURPOSE Set derivatives needed to determine a monotone piecewise 5C cubic Hermite interpolant to given data. Boundary values 6C are provided which are compatible with monotonicity. The 7C interpolant will have an extremum at each point where mono- 8C tonicity switches direction. (See DPCHIC if user control 9C is desired over boundary or switch conditions.) 10C***LIBRARY SLATEC (PCHIP) 11C***CATEGORY E1A 12C***TYPE DOUBLE PRECISION (PCHIM-S, DPCHIM-D) 13C***KEYWORDS CUBIC HERMITE INTERPOLATION, MONOTONE INTERPOLATION, 14C PCHIP, PIECEWISE CUBIC INTERPOLATION 15C***AUTHOR Fritsch, F. N., (LLNL) 16C Lawrence Livermore National Laboratory 17C P.O. Box 808 (L-316) 18C Livermore, CA 94550 19C FTS 532-4275, (510) 422-4275 20C***DESCRIPTION 21C 22C DPCHIM: Piecewise Cubic Hermite Interpolation to 23C Monotone data. 24C 25C Sets derivatives needed to determine a monotone piecewise cubic 26C Hermite interpolant to the data given in X and F. 27C 28C Default boundary conditions are provided which are compatible 29C with monotonicity. (See DPCHIC if user control of boundary con- 30C ditions is desired.) 31C 32C If the data are only piecewise monotonic, the interpolant will 33C have an extremum at each point where monotonicity switches direc- 34C tion. (See DPCHIC if user control is desired in such cases.) 35C 36C To facilitate two-dimensional applications, includes an increment 37C between successive values of the F- and D-arrays. 38C 39C The resulting piecewise cubic Hermite function may be evaluated 40C by DPCHFE or DPCHFD. 41C 42C ---------------------------------------------------------------------- 43C 44C Calling sequence: 45C 46C PARAMETER (INCFD = ...) 47C INTEGER N, IERR 48C DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N) 49C 50C CALL DPCHIM (N, X, F, D, INCFD, IERR) 51C 52C Parameters: 53C 54C N -- (input) number of data points. (Error return if N.LT.2 .) 55C If N=2, simply does linear interpolation. 56C 57C X -- (input) real*8 array of independent variable values. The 58C elements of X must be strictly increasing: 59C X(I-1) .LT. X(I), I = 2(1)N. 60C (Error return if not.) 61C 62C F -- (input) real*8 array of dependent variable values to be 63C interpolated. F(1+(I-1)*INCFD) is value corresponding to 64C X(I). DPCHIM is designed for monotonic data, but it will 65C work for any F-array. It will force extrema at points where 66C monotonicity switches direction. If some other treatment of 67C switch points is desired, DPCHIC should be used instead. 68C ----- 69C D -- (output) real*8 array of derivative values at the data 70C points. If the data are monotonic, these values will 71C determine a monotone cubic Hermite function. 72C The value corresponding to X(I) is stored in 73C D(1+(I-1)*INCFD), I=1(1)N. 74C No other entries in D are changed. 75C 76C INCFD -- (input) increment between successive values in F and D. 77C This argument is provided primarily for 2-D applications. 78C (Error return if INCFD.LT.1 .) 79C 80C IERR -- (output) error flag. 81C Normal return: 82C IERR = 0 (no errors). 83C Warning error: 84C IERR.GT.0 means that IERR switches in the direction 85C of monotonicity were detected. 86C "Recoverable" errors: 87C IERR = -1 if N.LT.2 . 88C IERR = -2 if INCFD.LT.1 . 89C IERR = -3 if the X-array is not strictly increasing. 90C (The D-array has not been changed in any of these cases.) 91C NOTE: The above errors are checked in the order listed, 92C and following arguments have **NOT** been validated. 93C 94C***REFERENCES 1. F. N. Fritsch and J. Butland, A method for construc- 95C ting local monotone piecewise cubic interpolants, SIAM 96C Journal on Scientific and Statistical Computing 5, 2 97C (June 1984), pp. 300-304. 98C 2. F. N. Fritsch and R. E. Carlson, Monotone piecewise 99C cubic interpolation, SIAM Journal on Numerical Ana- 100C lysis 17, 2 (April 1980), pp. 238-246. 101C***ROUTINES CALLED DPCHST, XERMSG 102C***REVISION HISTORY (YYMMDD) 103C 811103 DATE WRITTEN 104C 820201 1. Introduced DPCHST to reduce possible over/under- 105C flow problems. 106C 2. Rearranged derivative formula for same reason. 107C 820602 1. Modified end conditions to be continuous functions 108C of data when monotonicity switches in next interval. 109C 2. Modified formulas so end conditions are less prone 110C of over/underflow problems. 111C 820803 Minor cosmetic changes for release 1. 112C 870707 Corrected XERROR calls for d.p. name(s). 113C 870813 Updated Reference 1. 114C 890206 Corrected XERROR calls. 115C 890411 Added SAVE statements (Vers. 3.2). 116C 890531 Changed all specific intrinsics to generic. (WRB) 117C 890703 Corrected category record. (WRB) 118C 890831 Modified array declarations. (WRB) 119C 891006 Cosmetic changes to prologue. (WRB) 120C 891006 REVISION DATE from Version 3.2 121C 891214 Prologue converted to Version 4.0 format. (BAB) 122C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 123C 920429 Revised format and order of references. (WRB,FNF) 124C***END PROLOGUE DPCHIM 125C Programming notes: 126C 127C 1. The function DPCHST(ARG1,ARG2) is assumed to return zero if 128C either argument is zero, +1 if they are of the same sign, and 129C -1 if they are of opposite sign. 130C 2. To produce a single precision version, simply: 131C a. Change DPCHIM to PCHIM wherever it occurs, 132C b. Change DPCHST to PCHST wherever it occurs, 133C c. Change all references to the Fortran intrinsics to their 134C single precision equivalents, 135C d. Change the double precision declarations to real, and 136C e. Change the constants ZERO and THREE to single precision. 137C 138C DECLARE ARGUMENTS. 139C 140 INTEGER N, INCFD, IERR 141 DOUBLE PRECISION X(*), F(INCFD,*), D(INCFD,*) 142C 143C DECLARE LOCAL VARIABLES. 144C 145 INTEGER I, NLESS1 146 DOUBLE PRECISION DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE, 147 * H1, H2, HSUM, HSUMT3, THREE, W1, W2, ZERO 148 SAVE ZERO, THREE 149 DOUBLE PRECISION DPCHST 150 DATA ZERO /0.D0/, THREE/3.D0/ 151C 152C VALIDITY-CHECK ARGUMENTS. 153C 154C***FIRST EXECUTABLE STATEMENT DPCHIM 155 IF ( N.LT.2 ) GO TO 5001 156 IF ( INCFD.LT.1 ) GO TO 5002 157 DO 1 I = 2, N 158 IF ( X(I).LE.X(I-1) ) GO TO 5003 159 1 CONTINUE 160C 161C FUNCTION DEFINITION IS OK, GO ON. 162C 163 IERR = 0 164 NLESS1 = N - 1 165 H1 = X(2) - X(1) 166 DEL1 = (F(1,2) - F(1,1))/H1 167 DSAVE = DEL1 168C 169C SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. 170C 171 IF (NLESS1 .GT. 1) GO TO 10 172 D(1,1) = DEL1 173 D(1,N) = DEL1 174 GO TO 5000 175C 176C NORMAL CASE (N .GE. 3). 177C 178 10 CONTINUE 179 H2 = X(3) - X(2) 180 DEL2 = (F(1,3) - F(1,2))/H2 181C 182C SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE 183C SHAPE-PRESERVING. 184C 185 HSUM = H1 + H2 186 W1 = (H1 + HSUM)/HSUM 187 W2 = -H1/HSUM 188 D(1,1) = W1*DEL1 + W2*DEL2 189 IF ( DPCHST(D(1,1),DEL1) .LE. ZERO) THEN 190 D(1,1) = ZERO 191 ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO) THEN 192C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. 193 DMAX = THREE*DEL1 194 IF (ABS(D(1,1)) .GT. ABS(DMAX)) D(1,1) = DMAX 195 ENDIF 196C 197C LOOP THROUGH INTERIOR POINTS. 198C 199 DO 50 I = 2, NLESS1 200 IF (I .EQ. 2) GO TO 40 201C 202 H1 = H2 203 H2 = X(I+1) - X(I) 204 HSUM = H1 + H2 205 DEL1 = DEL2 206 DEL2 = (F(1,I+1) - F(1,I))/H2 207 40 CONTINUE 208C 209C SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. 210C 211 D(1,I) = ZERO 212 IF ( DPCHST(DEL1,DEL2) .LT. 0.) GO TO 42 213 IF ( DPCHST(DEL1,DEL2) .EQ. 0.) GO TO 41 214 GO TO 45 215C 216C COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY. 217C 218 41 CONTINUE 219 IF (DEL2 .EQ. ZERO) GO TO 50 220 IF ( DPCHST(DSAVE,DEL2) .LT. ZERO) IERR = IERR + 1 221 DSAVE = DEL2 222 GO TO 50 223C 224 42 CONTINUE 225 IERR = IERR + 1 226 DSAVE = DEL2 227 GO TO 50 228C 229C USE BRODLIE MODIFICATION OF BUTLAND FORMULA. 230C 231 45 CONTINUE 232 HSUMT3 = HSUM+HSUM+HSUM 233 W1 = (HSUM + H1)/HSUMT3 234 W2 = (HSUM + H2)/HSUMT3 235 DMAX = MAX( ABS(DEL1), ABS(DEL2) ) 236 DMIN = MIN( ABS(DEL1), ABS(DEL2) ) 237 DRAT1 = DEL1/DMAX 238 DRAT2 = DEL2/DMAX 239 D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2) 240C 241 50 CONTINUE 242C 243C SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE 244C SHAPE-PRESERVING. 245C 246 W1 = -H2/HSUM 247 W2 = (H2 + HSUM)/HSUM 248 D(1,N) = W1*DEL1 + W2*DEL2 249 IF ( DPCHST(D(1,N),DEL2) .LE. ZERO) THEN 250 D(1,N) = ZERO 251 ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO) THEN 252C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. 253 DMAX = THREE*DEL2 254 IF (ABS(D(1,N)) .GT. ABS(DMAX)) D(1,N) = DMAX 255 ENDIF 256C 257C NORMAL RETURN. 258C 259 5000 CONTINUE 260 RETURN 261C 262C ERROR RETURNS. 263C 264 5001 CONTINUE 265C N.LT.2 RETURN. 266 IERR = -1 267 CALL XERMSG ('SLATEC', 'DPCHIM', 268 + 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1) 269 RETURN 270C 271 5002 CONTINUE 272C INCFD.LT.1 RETURN. 273 IERR = -2 274 CALL XERMSG ('SLATEC', 'DPCHIM', 'INCREMENT LESS THAN ONE', IERR, 275 + 1) 276 RETURN 277C 278 5003 CONTINUE 279C X-ARRAY NOT STRICTLY INCREASING. 280 IERR = -3 281 CALL XERMSG ('SLATEC', 'DPCHIM', 282 + 'X-ARRAY NOT STRICTLY INCREASING', IERR, 1) 283 RETURN 284C------------- LAST LINE OF DPCHIM FOLLOWS ----------------------------- 285 END 286