1*DECK PCHIM 2 SUBROUTINE PCHIM (N, X, F, D, INCFD, IERR) 3C***BEGIN PROLOGUE PCHIM 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 PCHIC if user control is 9C desired over boundary or switch conditions.) 10C***LIBRARY SLATEC (PCHIP) 11C***CATEGORY E1A 12C***TYPE SINGLE 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 PCHIM: 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 PCHIC 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 PCHIC 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 PCHFE or PCHFD. 41C 42C ---------------------------------------------------------------------- 43C 44C Calling sequence: 45C 46C PARAMETER (INCFD = ...) 47C INTEGER N, IERR 48C REAL X(N), F(INCFD,N), D(INCFD,N) 49C 50C CALL PCHIM (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 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 array of dependent variable values to be inter- 63C polated. F(1+(I-1)*INCFD) is value corresponding to X(I). 64C PCHIM is designed for monotonic data, but it will work for 65C any F-array. It will force extrema at points where mono- 66C tonicity switches direction. If some other treatment of 67C switch points is desired, PCHIC should be used instead. 68C ----- 69C D -- (output) real array of derivative values at the data points. 70C If the data are monotonic, these values will determine a 71C 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 PCHST, XERMSG 102C***REVISION HISTORY (YYMMDD) 103C 811103 DATE WRITTEN 104C 820201 1. Introduced PCHST 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 870813 Updated Reference 1. 113C 890411 Added SAVE statements (Vers. 3.2). 114C 890531 Changed all specific intrinsics to generic. (WRB) 115C 890703 Corrected category record. (WRB) 116C 890831 Modified array declarations. (WRB) 117C 890831 REVISION DATE from Version 3.2 118C 891214 Prologue converted to Version 4.0 format. (BAB) 119C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 120C 920429 Revised format and order of references. (WRB,FNF) 121C***END PROLOGUE PCHIM 122C Programming notes: 123C 124C 1. The function PCHST(ARG1,ARG2) is assumed to return zero if 125C either argument is zero, +1 if they are of the same sign, and 126C -1 if they are of opposite sign. 127C 2. To produce a double precision version, simply: 128C a. Change PCHIM to DPCHIM wherever it occurs, 129C b. Change PCHST to DPCHST wherever it occurs, 130C c. Change all references to the Fortran intrinsics to their 131C double precision equivalents, 132C d. Change the real declarations to double precision, and 133C e. Change the constants ZERO and THREE to double precision. 134C 135C DECLARE ARGUMENTS. 136C 137 INTEGER N, INCFD, IERR 138 REAL X(*), F(INCFD,*), D(INCFD,*) 139C 140C DECLARE LOCAL VARIABLES. 141C 142 INTEGER I, NLESS1 143 REAL DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE, 144 * H1, H2, HSUM, HSUMT3, THREE, W1, W2, ZERO 145 SAVE ZERO, THREE 146 REAL PCHST 147 DATA ZERO /0./, THREE /3./ 148C 149C VALIDITY-CHECK ARGUMENTS. 150C 151C***FIRST EXECUTABLE STATEMENT PCHIM 152 IF ( N.LT.2 ) GO TO 5001 153 IF ( INCFD.LT.1 ) GO TO 5002 154 DO 1 I = 2, N 155 IF ( X(I).LE.X(I-1) ) GO TO 5003 156 1 CONTINUE 157C 158C FUNCTION DEFINITION IS OK, GO ON. 159C 160 IERR = 0 161 NLESS1 = N - 1 162 H1 = X(2) - X(1) 163 DEL1 = (F(1,2) - F(1,1))/H1 164 DSAVE = DEL1 165C 166C SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. 167C 168 IF (NLESS1 .GT. 1) GO TO 10 169 D(1,1) = DEL1 170 D(1,N) = DEL1 171 GO TO 5000 172C 173C NORMAL CASE (N .GE. 3). 174C 175 10 CONTINUE 176 H2 = X(3) - X(2) 177 DEL2 = (F(1,3) - F(1,2))/H2 178C 179C SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE 180C SHAPE-PRESERVING. 181C 182 HSUM = H1 + H2 183 W1 = (H1 + HSUM)/HSUM 184 W2 = -H1/HSUM 185 D(1,1) = W1*DEL1 + W2*DEL2 186 IF ( PCHST(D(1,1),DEL1) .LE. ZERO) THEN 187 D(1,1) = ZERO 188 ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO) THEN 189C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. 190 DMAX = THREE*DEL1 191 IF (ABS(D(1,1)) .GT. ABS(DMAX)) D(1,1) = DMAX 192 ENDIF 193C 194C LOOP THROUGH INTERIOR POINTS. 195C 196 DO 50 I = 2, NLESS1 197 IF (I .EQ. 2) GO TO 40 198C 199 H1 = H2 200 H2 = X(I+1) - X(I) 201 HSUM = H1 + H2 202 DEL1 = DEL2 203 DEL2 = (F(1,I+1) - F(1,I))/H2 204 40 CONTINUE 205C 206C SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. 207C 208 D(1,I) = ZERO 209 IF ( PCHST(DEL1,DEL2) ) 42, 41, 45 210C 211C COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY. 212C 213 41 CONTINUE 214 IF (DEL2 .EQ. ZERO) GO TO 50 215 IF ( PCHST(DSAVE,DEL2) .LT. ZERO) IERR = IERR + 1 216 DSAVE = DEL2 217 GO TO 50 218C 219 42 CONTINUE 220 IERR = IERR + 1 221 DSAVE = DEL2 222 GO TO 50 223C 224C USE BRODLIE MODIFICATION OF BUTLAND FORMULA. 225C 226 45 CONTINUE 227 HSUMT3 = HSUM+HSUM+HSUM 228 W1 = (HSUM + H1)/HSUMT3 229 W2 = (HSUM + H2)/HSUMT3 230 DMAX = MAX( ABS(DEL1), ABS(DEL2) ) 231 DMIN = MIN( ABS(DEL1), ABS(DEL2) ) 232 DRAT1 = DEL1/DMAX 233 DRAT2 = DEL2/DMAX 234 D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2) 235C 236 50 CONTINUE 237C 238C SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE 239C SHAPE-PRESERVING. 240C 241 W1 = -H2/HSUM 242 W2 = (H2 + HSUM)/HSUM 243 D(1,N) = W1*DEL1 + W2*DEL2 244 IF ( PCHST(D(1,N),DEL2) .LE. ZERO) THEN 245 D(1,N) = ZERO 246 ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO) THEN 247C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. 248 DMAX = THREE*DEL2 249 IF (ABS(D(1,N)) .GT. ABS(DMAX)) D(1,N) = DMAX 250 ENDIF 251C 252C NORMAL RETURN. 253C 254 5000 CONTINUE 255 RETURN 256C 257C ERROR RETURNS. 258C 259 5001 CONTINUE 260C N.LT.2 RETURN. 261 IERR = -1 262 CALL XERMSG ('SLATEC', 'PCHIM', 263 + 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1) 264 RETURN 265C 266 5002 CONTINUE 267C INCFD.LT.1 RETURN. 268 IERR = -2 269 CALL XERMSG ('SLATEC', 'PCHIM', 'INCREMENT LESS THAN ONE', IERR, 270 + 1) 271 RETURN 272C 273 5003 CONTINUE 274C X-ARRAY NOT STRICTLY INCREASING. 275 IERR = -3 276 CALL XERMSG ('SLATEC', 'PCHIM', 'X-ARRAY NOT STRICTLY INCREASING' 277 + , IERR, 1) 278 RETURN 279C------------- LAST LINE OF PCHIM FOLLOWS ------------------------------ 280 END 281