00001 *DECK DPCHIM 00002 SUBROUTINE DPCHIM (N, X, F, D, INCFD, IERR) 00003 C***BEGIN PROLOGUE DPCHIM 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 DPCHIC if user control 00009 C is desired over boundary or switch conditions.) 00010 C***LIBRARY SLATEC (PCHIP) 00011 C***CATEGORY E1A 00012 C***TYPE DOUBLE 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 DPCHIM: 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 DPCHIC 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 DPCHIC 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 DPCHFE or DPCHFD. 00041 C 00042 C ---------------------------------------------------------------------- 00043 C 00044 C Calling sequence: 00045 C 00046 C PARAMETER (INCFD = ...) 00047 C INTEGER N, IERR 00048 C DOUBLE PRECISION X(N), F(INCFD,N), D(INCFD,N) 00049 C 00050 C CALL DPCHIM (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*8 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*8 array of dependent variable values to be 00063 C interpolated. F(1+(I-1)*INCFD) is value corresponding to 00064 C X(I). DPCHIM is designed for monotonic data, but it will 00065 C work for any F-array. It will force extrema at points where 00066 C monotonicity switches direction. If some other treatment of 00067 C switch points is desired, DPCHIC should be used instead. 00068 C ----- 00069 C D -- (output) real*8 array of derivative values at the data 00070 C points. If the data are monotonic, these values will 00071 C determine 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 DPCHST, XERMSG 00102 C***REVISION HISTORY (YYMMDD) 00103 C 811103 DATE WRITTEN 00104 C 820201 1. Introduced DPCHST 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 870707 Corrected XERROR calls for d.p. name(s). 00113 C 870813 Updated Reference 1. 00114 C 890206 Corrected XERROR calls. 00115 C 890411 Added SAVE statements (Vers. 3.2). 00116 C 890531 Changed all specific intrinsics to generic. (WRB) 00117 C 890703 Corrected category record. (WRB) 00118 C 890831 Modified array declarations. (WRB) 00119 C 891006 Cosmetic changes to prologue. (WRB) 00120 C 891006 REVISION DATE from Version 3.2 00121 C 891214 Prologue converted to Version 4.0 format. (BAB) 00122 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 00123 C 920429 Revised format and order of references. (WRB,FNF) 00124 C***END PROLOGUE DPCHIM 00125 C Programming notes: 00126 C 00127 C 1. The function DPCHST(ARG1,ARG2) is assumed to return zero if 00128 C either argument is zero, +1 if they are of the same sign, and 00129 C -1 if they are of opposite sign. 00130 C 2. To produce a single precision version, simply: 00131 C a. Change DPCHIM to PCHIM wherever it occurs, 00132 C b. Change DPCHST to PCHST wherever it occurs, 00133 C c. Change all references to the Fortran intrinsics to their 00134 C single precision equivalents, 00135 C d. Change the double precision declarations to real, and 00136 C e. Change the constants ZERO and THREE to single precision. 00137 C 00138 C DECLARE ARGUMENTS. 00139 C 00140 INTEGER N, INCFD, IERR 00141 DOUBLE PRECISION X(*), F(INCFD,*), D(INCFD,*) 00142 C 00143 C DECLARE LOCAL VARIABLES. 00144 C 00145 INTEGER I, NLESS1 00146 DOUBLE PRECISION DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE, 00147 * H1, H2, HSUM, HSUMT3, THREE, W1, W2, ZERO 00148 SAVE ZERO, THREE 00149 DOUBLE PRECISION DPCHST 00150 DATA ZERO /0.D0/, THREE/3.D0/ 00151 C 00152 C VALIDITY-CHECK ARGUMENTS. 00153 C 00154 C***FIRST EXECUTABLE STATEMENT DPCHIM 00155 IF ( N.LT.2 ) GO TO 5001 00156 IF ( INCFD.LT.1 ) GO TO 5002 00157 DO 1 I = 2, N 00158 IF ( X(I).LE.X(I-1) ) GO TO 5003 00159 1 CONTINUE 00160 C 00161 C FUNCTION DEFINITION IS OK, GO ON. 00162 C 00163 IERR = 0 00164 NLESS1 = N - 1 00165 H1 = X(2) - X(1) 00166 DEL1 = (F(1,2) - F(1,1))/H1 00167 DSAVE = DEL1 00168 C 00169 C SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. 00170 C 00171 IF (NLESS1 .GT. 1) GO TO 10 00172 D(1,1) = DEL1 00173 D(1,N) = DEL1 00174 GO TO 5000 00175 C 00176 C NORMAL CASE (N .GE. 3). 00177 C 00178 10 CONTINUE 00179 H2 = X(3) - X(2) 00180 DEL2 = (F(1,3) - F(1,2))/H2 00181 C 00182 C SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE 00183 C SHAPE-PRESERVING. 00184 C 00185 HSUM = H1 + H2 00186 W1 = (H1 + HSUM)/HSUM 00187 W2 = -H1/HSUM 00188 D(1,1) = W1*DEL1 + W2*DEL2 00189 IF ( DPCHST(D(1,1),DEL1) .LE. ZERO) THEN 00190 D(1,1) = ZERO 00191 ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO) THEN 00192 C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. 00193 DMAX = THREE*DEL1 00194 IF (ABS(D(1,1)) .GT. ABS(DMAX)) D(1,1) = DMAX 00195 ENDIF 00196 C 00197 C LOOP THROUGH INTERIOR POINTS. 00198 C 00199 DO 50 I = 2, NLESS1 00200 IF (I .EQ. 2) GO TO 40 00201 C 00202 H1 = H2 00203 H2 = X(I+1) - X(I) 00204 HSUM = H1 + H2 00205 DEL1 = DEL2 00206 DEL2 = (F(1,I+1) - F(1,I))/H2 00207 40 CONTINUE 00208 C 00209 C SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. 00210 C 00211 D(1,I) = ZERO 00212 IF ( DPCHST(DEL1,DEL2) .LT. 0.) GO TO 42 00213 IF ( DPCHST(DEL1,DEL2) .EQ. 0.) GO TO 41 00214 GO TO 45 00215 C 00216 C COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY. 00217 C 00218 41 CONTINUE 00219 IF (DEL2 .EQ. ZERO) GO TO 50 00220 IF ( DPCHST(DSAVE,DEL2) .LT. ZERO) IERR = IERR + 1 00221 DSAVE = DEL2 00222 GO TO 50 00223 C 00224 42 CONTINUE 00225 IERR = IERR + 1 00226 DSAVE = DEL2 00227 GO TO 50 00228 C 00229 C USE BRODLIE MODIFICATION OF BUTLAND FORMULA. 00230 C 00231 45 CONTINUE 00232 HSUMT3 = HSUM+HSUM+HSUM 00233 W1 = (HSUM + H1)/HSUMT3 00234 W2 = (HSUM + H2)/HSUMT3 00235 DMAX = MAX( ABS(DEL1), ABS(DEL2) ) 00236 DMIN = MIN( ABS(DEL1), ABS(DEL2) ) 00237 DRAT1 = DEL1/DMAX 00238 DRAT2 = DEL2/DMAX 00239 D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2) 00240 C 00241 50 CONTINUE 00242 C 00243 C SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE 00244 C SHAPE-PRESERVING. 00245 C 00246 W1 = -H2/HSUM 00247 W2 = (H2 + HSUM)/HSUM 00248 D(1,N) = W1*DEL1 + W2*DEL2 00249 IF ( DPCHST(D(1,N),DEL2) .LE. ZERO) THEN 00250 D(1,N) = ZERO 00251 ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO) THEN 00252 C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. 00253 DMAX = THREE*DEL2 00254 IF (ABS(D(1,N)) .GT. ABS(DMAX)) D(1,N) = DMAX 00255 ENDIF 00256 C 00257 C NORMAL RETURN. 00258 C 00259 5000 CONTINUE 00260 RETURN 00261 C 00262 C ERROR RETURNS. 00263 C 00264 5001 CONTINUE 00265 C N.LT.2 RETURN. 00266 IERR = -1 00267 CALL XERMSG ('SLATEC', 'DPCHIM', 00268 + 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1) 00269 RETURN 00270 C 00271 5002 CONTINUE 00272 C INCFD.LT.1 RETURN. 00273 IERR = -2 00274 CALL XERMSG ('SLATEC', 'DPCHIM', 'INCREMENT LESS THAN ONE', IERR, 00275 + 1) 00276 RETURN 00277 C 00278 5003 CONTINUE 00279 C X-ARRAY NOT STRICTLY INCREASING. 00280 IERR = -3 00281 CALL XERMSG ('SLATEC', 'DPCHIM', 00282 + 'X-ARRAY NOT STRICTLY INCREASING', IERR, 1) 00283 RETURN 00284 C------------- LAST LINE OF DPCHIM FOLLOWS ----------------------------- 00285 END
1.7.1
