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