1*DECK DPCHCI 2 SUBROUTINE DPCHCI (N, H, SLOPE, D, INCFD) 3C***BEGIN PROLOGUE DPCHCI 4C***SUBSIDIARY 5C***PURPOSE Set interior derivatives for DPCHIC 6C***LIBRARY SLATEC (PCHIP) 7C***TYPE DOUBLE PRECISION (PCHCI-S, DPCHCI-D) 8C***AUTHOR Fritsch, F. N., (LLNL) 9C***DESCRIPTION 10C 11C DPCHCI: DPCHIC Initial Derivative Setter. 12C 13C Called by DPCHIC to set derivatives needed to determine a monotone 14C piecewise cubic Hermite interpolant to the data. 15C 16C Default boundary conditions are provided which are compatible 17C with monotonicity. If the data are only piecewise monotonic, the 18C interpolant will have an extremum at each point where monotonicity 19C switches direction. 20C 21C To facilitate two-dimensional applications, includes an increment 22C between successive values of the D-array. 23C 24C The resulting piecewise cubic Hermite function should be identical 25C (within roundoff error) to that produced by DPCHIM. 26C 27C ---------------------------------------------------------------------- 28C 29C Calling sequence: 30C 31C PARAMETER (INCFD = ...) 32C INTEGER N 33C DOUBLE PRECISION H(N), SLOPE(N), D(INCFD,N) 34C 35C CALL DPCHCI (N, H, SLOPE, D, INCFD) 36C 37C Parameters: 38C 39C N -- (input) number of data points. 40C If N=2, simply does linear interpolation. 41C 42C H -- (input) real*8 array of interval lengths. 43C SLOPE -- (input) real*8 array of data slopes. 44C If the data are (X(I),Y(I)), I=1(1)N, then these inputs are: 45C H(I) = X(I+1)-X(I), 46C SLOPE(I) = (Y(I+1)-Y(I))/H(I), I=1(1)N-1. 47C 48C D -- (output) real*8 array of derivative values at data points. 49C If the data are monotonic, these values will determine a 50C a monotone cubic Hermite function. 51C The value corresponding to X(I) is stored in 52C D(1+(I-1)*INCFD), I=1(1)N. 53C No other entries in D are changed. 54C 55C INCFD -- (input) increment between successive values in D. 56C This argument is provided primarily for 2-D applications. 57C 58C ------- 59C WARNING: This routine does no validity-checking of arguments. 60C ------- 61C 62C Fortran intrinsics used: ABS, MAX, MIN. 63C 64C***SEE ALSO DPCHIC 65C***ROUTINES CALLED DPCHST 66C***REVISION HISTORY (YYMMDD) 67C 820218 DATE WRITTEN 68C 820601 Modified end conditions to be continuous functions of 69C data when monotonicity switches in next interval. 70C 820602 1. Modified formulas so end conditions are less prone 71C to over/underflow problems. 72C 2. Minor modification to HSUM calculation. 73C 820805 Converted to SLATEC library version. 74C 870813 Minor cosmetic changes. 75C 890411 Added SAVE statements (Vers. 3.2). 76C 890531 Changed all specific intrinsics to generic. (WRB) 77C 890831 Modified array declarations. (WRB) 78C 890831 REVISION DATE from Version 3.2 79C 891214 Prologue converted to Version 4.0 format. (BAB) 80C 900328 Added TYPE section. (WRB) 81C 910408 Updated AUTHOR section in prologue. (WRB) 82C 930503 Improved purpose. (FNF) 83C***END PROLOGUE DPCHCI 84C 85C Programming notes: 86C 1. The function DPCHST(ARG1,ARG2) is assumed to return zero if 87C either argument is zero, +1 if they are of the same sign, and 88C -1 if they are of opposite sign. 89C**End 90C 91C DECLARE ARGUMENTS. 92C 93 INTEGER N, INCFD 94 DOUBLE PRECISION H(*), SLOPE(*), D(INCFD,*) 95C 96C DECLARE LOCAL VARIABLES. 97C 98 INTEGER I, NLESS1 99 DOUBLE PRECISION DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, HSUM, 100 * HSUMT3, THREE, W1, W2, ZERO 101 SAVE ZERO, THREE 102 DOUBLE PRECISION DPCHST 103C 104C INITIALIZE. 105C 106 DATA ZERO /0.D0/, THREE/3.D0/ 107C***FIRST EXECUTABLE STATEMENT DPCHCI 108 NLESS1 = N - 1 109 DEL1 = SLOPE(1) 110C 111C SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. 112C 113 IF (NLESS1 .GT. 1) GO TO 10 114 D(1,1) = DEL1 115 D(1,N) = DEL1 116 GO TO 5000 117C 118C NORMAL CASE (N .GE. 3). 119C 120 10 CONTINUE 121 DEL2 = SLOPE(2) 122C 123C SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE 124C SHAPE-PRESERVING. 125C 126 HSUM = H(1) + H(2) 127 W1 = (H(1) + HSUM)/HSUM 128 W2 = -H(1)/HSUM 129 D(1,1) = W1*DEL1 + W2*DEL2 130 IF ( DPCHST(D(1,1),DEL1) .LE. ZERO) THEN 131 D(1,1) = ZERO 132 ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO) THEN 133C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. 134 DMAX = THREE*DEL1 135 IF (ABS(D(1,1)) .GT. ABS(DMAX)) D(1,1) = DMAX 136 ENDIF 137C 138C LOOP THROUGH INTERIOR POINTS. 139C 140 DO 50 I = 2, NLESS1 141 IF (I .EQ. 2) GO TO 40 142C 143 HSUM = H(I-1) + H(I) 144 DEL1 = DEL2 145 DEL2 = SLOPE(I) 146 40 CONTINUE 147C 148C SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. 149C 150 D(1,I) = ZERO 151 IF ( DPCHST(DEL1,DEL2) .LE. ZERO) GO TO 50 152C 153C USE BRODLIE MODIFICATION OF BUTLAND FORMULA. 154C 155 HSUMT3 = HSUM+HSUM+HSUM 156 W1 = (HSUM + H(I-1))/HSUMT3 157 W2 = (HSUM + H(I) )/HSUMT3 158 DMAX = MAX( ABS(DEL1), ABS(DEL2) ) 159 DMIN = MIN( ABS(DEL1), ABS(DEL2) ) 160 DRAT1 = DEL1/DMAX 161 DRAT2 = DEL2/DMAX 162 D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2) 163C 164 50 CONTINUE 165C 166C SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE 167C SHAPE-PRESERVING. 168C 169 W1 = -H(N-1)/HSUM 170 W2 = (H(N-1) + HSUM)/HSUM 171 D(1,N) = W1*DEL1 + W2*DEL2 172 IF ( DPCHST(D(1,N),DEL2) .LE. ZERO) THEN 173 D(1,N) = ZERO 174 ELSE IF ( DPCHST(DEL1,DEL2) .LT. ZERO) THEN 175C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. 176 DMAX = THREE*DEL2 177 IF (ABS(D(1,N)) .GT. ABS(DMAX)) D(1,N) = DMAX 178 ENDIF 179C 180C NORMAL RETURN. 181C 182 5000 CONTINUE 183 RETURN 184C------------- LAST LINE OF DPCHCI FOLLOWS ----------------------------- 185 END 186