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