00001 *DECK PCHIM 00002 SUBROUTINE PCHIM (N, X, F, D, INCFD, IERR) 00003 C***BEGIN PROLOGUE PCHIM 00004 C***PURPOSE Set derivatives needed to determine a monotone piecewise 00005 C cubic Hermite interpolant to given data. Boundary values 00006 C are provided which are compatible with monotonicity. The 00007 C interpolant will have an extremum at each point where mono- 00008 C tonicity switches direction. (See PCHIC if user control is 00009 C desired over boundary or switch conditions.) 00010 C***LIBRARY SLATEC (PCHIP) 00011 C***CATEGORY E1A 00012 C***TYPE SINGLE PRECISION (PCHIM-S, DPCHIM-D) 00013 C***KEYWORDS CUBIC HERMITE INTERPOLATION, MONOTONE INTERPOLATION, 00014 C PCHIP, PIECEWISE CUBIC INTERPOLATION 00015 C***AUTHOR Fritsch, F. N., (LLNL) 00016 C Lawrence Livermore National Laboratory 00017 C P.O. Box 808 (L-316) 00018 C Livermore, CA 94550 00019 C FTS 532-4275, (510) 422-4275 00020 C***DESCRIPTION 00021 C 00022 C PCHIM: Piecewise Cubic Hermite Interpolation to 00023 C Monotone data. 00024 C 00025 C Sets derivatives needed to determine a monotone piecewise cubic 00026 C Hermite interpolant to the data given in X and F. 00027 C 00028 C Default boundary conditions are provided which are compatible 00029 C with monotonicity. (See PCHIC if user control of boundary con- 00030 C ditions is desired.) 00031 C 00032 C If the data are only piecewise monotonic, the interpolant will 00033 C have an extremum at each point where monotonicity switches direc- 00034 C tion. (See PCHIC if user control is desired in such cases.) 00035 C 00036 C To facilitate two-dimensional applications, includes an increment 00037 C between successive values of the F- and D-arrays. 00038 C 00039 C The resulting piecewise cubic Hermite function may be evaluated 00040 C by PCHFE or PCHFD. 00041 C 00042 C ---------------------------------------------------------------------- 00043 C 00044 C Calling sequence: 00045 C 00046 C PARAMETER (INCFD = ...) 00047 C INTEGER N, IERR 00048 C REAL X(N), F(INCFD,N), D(INCFD,N) 00049 C 00050 C CALL PCHIM (N, X, F, D, INCFD, IERR) 00051 C 00052 C Parameters: 00053 C 00054 C N -- (input) number of data points. (Error return if N.LT.2 .) 00055 C If N=2, simply does linear interpolation. 00056 C 00057 C X -- (input) real array of independent variable values. The 00058 C elements of X must be strictly increasing: 00059 C X(I-1) .LT. X(I), I = 2(1)N. 00060 C (Error return if not.) 00061 C 00062 C F -- (input) real array of dependent variable values to be inter- 00063 C polated. F(1+(I-1)*INCFD) is value corresponding to X(I). 00064 C PCHIM is designed for monotonic data, but it will work for 00065 C any F-array. It will force extrema at points where mono- 00066 C tonicity switches direction. If some other treatment of 00067 C switch points is desired, PCHIC should be used instead. 00068 C ----- 00069 C D -- (output) real array of derivative values at the data points. 00070 C If the data are monotonic, these values will determine a 00071 C a monotone cubic Hermite function. 00072 C The value corresponding to X(I) is stored in 00073 C D(1+(I-1)*INCFD), I=1(1)N. 00074 C No other entries in D are changed. 00075 C 00076 C INCFD -- (input) increment between successive values in F and D. 00077 C This argument is provided primarily for 2-D applications. 00078 C (Error return if INCFD.LT.1 .) 00079 C 00080 C IERR -- (output) error flag. 00081 C Normal return: 00082 C IERR = 0 (no errors). 00083 C Warning error: 00084 C IERR.GT.0 means that IERR switches in the direction 00085 C of monotonicity were detected. 00086 C "Recoverable" errors: 00087 C IERR = -1 if N.LT.2 . 00088 C IERR = -2 if INCFD.LT.1 . 00089 C IERR = -3 if the X-array is not strictly increasing. 00090 C (The D-array has not been changed in any of these cases.) 00091 C NOTE: The above errors are checked in the order listed, 00092 C and following arguments have **NOT** been validated. 00093 C 00094 C***REFERENCES 1. F. N. Fritsch and J. Butland, A method for construc- 00095 C ting local monotone piecewise cubic interpolants, SIAM 00096 C Journal on Scientific and Statistical Computing 5, 2 00097 C (June 1984), pp. 300-304. 00098 C 2. F. N. Fritsch and R. E. Carlson, Monotone piecewise 00099 C cubic interpolation, SIAM Journal on Numerical Ana- 00100 C lysis 17, 2 (April 1980), pp. 238-246. 00101 C***ROUTINES CALLED PCHST, XERMSG 00102 C***REVISION HISTORY (YYMMDD) 00103 C 811103 DATE WRITTEN 00104 C 820201 1. Introduced PCHST to reduce possible over/under- 00105 C flow problems. 00106 C 2. Rearranged derivative formula for same reason. 00107 C 820602 1. Modified end conditions to be continuous functions 00108 C of data when monotonicity switches in next interval. 00109 C 2. Modified formulas so end conditions are less prone 00110 C of over/underflow problems. 00111 C 820803 Minor cosmetic changes for release 1. 00112 C 870813 Updated Reference 1. 00113 C 890411 Added SAVE statements (Vers. 3.2). 00114 C 890531 Changed all specific intrinsics to generic. (WRB) 00115 C 890703 Corrected category record. (WRB) 00116 C 890831 Modified array declarations. (WRB) 00117 C 890831 REVISION DATE from Version 3.2 00118 C 891214 Prologue converted to Version 4.0 format. (BAB) 00119 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 00120 C 920429 Revised format and order of references. (WRB,FNF) 00121 C***END PROLOGUE PCHIM 00122 C Programming notes: 00123 C 00124 C 1. The function PCHST(ARG1,ARG2) is assumed to return zero if 00125 C either argument is zero, +1 if they are of the same sign, and 00126 C -1 if they are of opposite sign. 00127 C 2. To produce a double precision version, simply: 00128 C a. Change PCHIM to DPCHIM wherever it occurs, 00129 C b. Change PCHST to DPCHST wherever it occurs, 00130 C c. Change all references to the Fortran intrinsics to their 00131 C double precision equivalents, 00132 C d. Change the real declarations to double precision, and 00133 C e. Change the constants ZERO and THREE to double precision. 00134 C 00135 C DECLARE ARGUMENTS. 00136 C 00137 INTEGER N, INCFD, IERR 00138 REAL X(*), F(INCFD,*), D(INCFD,*) 00139 C 00140 C DECLARE LOCAL VARIABLES. 00141 C 00142 INTEGER I, NLESS1 00143 REAL DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE, 00144 * H1, H2, HSUM, HSUMT3, THREE, W1, W2, ZERO 00145 SAVE ZERO, THREE 00146 REAL PCHST 00147 DATA ZERO /0./, THREE /3./ 00148 C 00149 C VALIDITY-CHECK ARGUMENTS. 00150 C 00151 C***FIRST EXECUTABLE STATEMENT PCHIM 00152 IF ( N.LT.2 ) GO TO 5001 00153 IF ( INCFD.LT.1 ) GO TO 5002 00154 DO 1 I = 2, N 00155 IF ( X(I).LE.X(I-1) ) GO TO 5003 00156 1 CONTINUE 00157 C 00158 C FUNCTION DEFINITION IS OK, GO ON. 00159 C 00160 IERR = 0 00161 NLESS1 = N - 1 00162 H1 = X(2) - X(1) 00163 DEL1 = (F(1,2) - F(1,1))/H1 00164 DSAVE = DEL1 00165 C 00166 C SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. 00167 C 00168 IF (NLESS1 .GT. 1) GO TO 10 00169 D(1,1) = DEL1 00170 D(1,N) = DEL1 00171 GO TO 5000 00172 C 00173 C NORMAL CASE (N .GE. 3). 00174 C 00175 10 CONTINUE 00176 H2 = X(3) - X(2) 00177 DEL2 = (F(1,3) - F(1,2))/H2 00178 C 00179 C SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE 00180 C SHAPE-PRESERVING. 00181 C 00182 HSUM = H1 + H2 00183 W1 = (H1 + HSUM)/HSUM 00184 W2 = -H1/HSUM 00185 D(1,1) = W1*DEL1 + W2*DEL2 00186 IF ( PCHST(D(1,1),DEL1) .LE. ZERO) THEN 00187 D(1,1) = ZERO 00188 ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO) THEN 00189 C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. 00190 DMAX = THREE*DEL1 00191 IF (ABS(D(1,1)) .GT. ABS(DMAX)) D(1,1) = DMAX 00192 ENDIF 00193 C 00194 C LOOP THROUGH INTERIOR POINTS. 00195 C 00196 DO 50 I = 2, NLESS1 00197 IF (I .EQ. 2) GO TO 40 00198 C 00199 H1 = H2 00200 H2 = X(I+1) - X(I) 00201 HSUM = H1 + H2 00202 DEL1 = DEL2 00203 DEL2 = (F(1,I+1) - F(1,I))/H2 00204 40 CONTINUE 00205 C 00206 C SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. 00207 C 00208 D(1,I) = ZERO 00209 IF ( PCHST(DEL1,DEL2) ) 42, 41, 45 00210 C 00211 C COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY. 00212 C 00213 41 CONTINUE 00214 IF (DEL2 .EQ. ZERO) GO TO 50 00215 IF ( PCHST(DSAVE,DEL2) .LT. ZERO) IERR = IERR + 1 00216 DSAVE = DEL2 00217 GO TO 50 00218 C 00219 42 CONTINUE 00220 IERR = IERR + 1 00221 DSAVE = DEL2 00222 GO TO 50 00223 C 00224 C USE BRODLIE MODIFICATION OF BUTLAND FORMULA. 00225 C 00226 45 CONTINUE 00227 HSUMT3 = HSUM+HSUM+HSUM 00228 W1 = (HSUM + H1)/HSUMT3 00229 W2 = (HSUM + H2)/HSUMT3 00230 DMAX = MAX( ABS(DEL1), ABS(DEL2) ) 00231 DMIN = MIN( ABS(DEL1), ABS(DEL2) ) 00232 DRAT1 = DEL1/DMAX 00233 DRAT2 = DEL2/DMAX 00234 D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2) 00235 C 00236 50 CONTINUE 00237 C 00238 C SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE 00239 C SHAPE-PRESERVING. 00240 C 00241 W1 = -H2/HSUM 00242 W2 = (H2 + HSUM)/HSUM 00243 D(1,N) = W1*DEL1 + W2*DEL2 00244 IF ( PCHST(D(1,N),DEL2) .LE. ZERO) THEN 00245 D(1,N) = ZERO 00246 ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO) THEN 00247 C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. 00248 DMAX = THREE*DEL2 00249 IF (ABS(D(1,N)) .GT. ABS(DMAX)) D(1,N) = DMAX 00250 ENDIF 00251 C 00252 C NORMAL RETURN. 00253 C 00254 5000 CONTINUE 00255 RETURN 00256 C 00257 C ERROR RETURNS. 00258 C 00259 5001 CONTINUE 00260 C N.LT.2 RETURN. 00261 IERR = -1 00262 CALL XERMSG ('SLATEC', 'PCHIM', 00263 + 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1) 00264 RETURN 00265 C 00266 5002 CONTINUE 00267 C INCFD.LT.1 RETURN. 00268 IERR = -2 00269 CALL XERMSG ('SLATEC', 'PCHIM', 'INCREMENT LESS THAN ONE', IERR, 00270 + 1) 00271 RETURN 00272 C 00273 5003 CONTINUE 00274 C X-ARRAY NOT STRICTLY INCREASING. 00275 IERR = -3 00276 CALL XERMSG ('SLATEC', 'PCHIM', 'X-ARRAY NOT STRICTLY INCREASING' 00277 + , IERR, 1) 00278 RETURN 00279 C------------- LAST LINE OF PCHIM FOLLOWS ------------------------------ 00280 END
1.7.1
