1*DECK PCHCI 2 SUBROUTINE PCHCI (N, H, SLOPE, D, INCFD) 3C***BEGIN PROLOGUE PCHCI 4C***SUBSIDIARY 5C***PURPOSE Set interior derivatives for PCHIC 6C***LIBRARY SLATEC (PCHIP) 7C***TYPE SINGLE PRECISION (PCHCI-S, DPCHCI-D) 8C***AUTHOR Fritsch, F. N., (LLNL) 9C***DESCRIPTION 10C 11C PCHCI: PCHIC Initial Derivative Setter. 12C 13C Called by PCHIC 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 PCHIM. 26C 27C ---------------------------------------------------------------------- 28C 29C Calling sequence: 30C 31C PARAMETER (INCFD = ...) 32C INTEGER N 33C REAL H(N), SLOPE(N), D(INCFD,N) 34C 35C CALL PCHCI (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 array of interval lengths. 43C SLOPE -- (input) real 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 array of derivative values at the 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 PCHIC 65C***ROUTINES CALLED PCHST 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 890411 Added SAVE statements (Vers. 3.2). 75C 890531 Changed all specific intrinsics to generic. (WRB) 76C 890831 Modified array declarations. (WRB) 77C 890831 REVISION DATE from Version 3.2 78C 891214 Prologue converted to Version 4.0 format. (BAB) 79C 900328 Added TYPE section. (WRB) 80C 910408 Updated AUTHOR section in prologue. (WRB) 81C 930503 Improved purpose. (FNF) 82C***END PROLOGUE PCHCI 83C 84C Programming notes: 85C 1. The function PCHST(ARG1,ARG2) is assumed to return zero if 86C either argument is zero, +1 if they are of the same sign, and 87C -1 if they are of opposite sign. 88C**End 89C 90C DECLARE ARGUMENTS. 91C 92 INTEGER N, INCFD 93 REAL H(*), SLOPE(*), D(INCFD,*) 94C 95C DECLARE LOCAL VARIABLES. 96C 97 INTEGER I, NLESS1 98 REAL DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, HSUM, HSUMT3, THREE, 99 * W1, W2, ZERO 100 SAVE ZERO, THREE 101 REAL PCHST 102C 103C INITIALIZE. 104C 105 DATA ZERO /0./, THREE /3./ 106C***FIRST EXECUTABLE STATEMENT PCHCI 107 NLESS1 = N - 1 108 DEL1 = SLOPE(1) 109C 110C SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. 111C 112 IF (NLESS1 .GT. 1) GO TO 10 113 D(1,1) = DEL1 114 D(1,N) = DEL1 115 GO TO 5000 116C 117C NORMAL CASE (N .GE. 3). 118C 119 10 CONTINUE 120 DEL2 = SLOPE(2) 121C 122C SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE 123C SHAPE-PRESERVING. 124C 125 HSUM = H(1) + H(2) 126 W1 = (H(1) + HSUM)/HSUM 127 W2 = -H(1)/HSUM 128 D(1,1) = W1*DEL1 + W2*DEL2 129 IF ( PCHST(D(1,1),DEL1) .LE. ZERO) THEN 130 D(1,1) = ZERO 131 ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO) THEN 132C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. 133 DMAX = THREE*DEL1 134 IF (ABS(D(1,1)) .GT. ABS(DMAX)) D(1,1) = DMAX 135 ENDIF 136C 137C LOOP THROUGH INTERIOR POINTS. 138C 139 DO 50 I = 2, NLESS1 140 IF (I .EQ. 2) GO TO 40 141C 142 HSUM = H(I-1) + H(I) 143 DEL1 = DEL2 144 DEL2 = SLOPE(I) 145 40 CONTINUE 146C 147C SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. 148C 149 D(1,I) = ZERO 150 IF ( PCHST(DEL1,DEL2) .LE. ZERO) GO TO 50 151C 152C USE BRODLIE MODIFICATION OF BUTLAND FORMULA. 153C 154 HSUMT3 = HSUM+HSUM+HSUM 155 W1 = (HSUM + H(I-1))/HSUMT3 156 W2 = (HSUM + H(I) )/HSUMT3 157 DMAX = MAX( ABS(DEL1), ABS(DEL2) ) 158 DMIN = MIN( ABS(DEL1), ABS(DEL2) ) 159 DRAT1 = DEL1/DMAX 160 DRAT2 = DEL2/DMAX 161 D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2) 162C 163 50 CONTINUE 164C 165C SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE 166C SHAPE-PRESERVING. 167C 168 W1 = -H(N-1)/HSUM 169 W2 = (H(N-1) + HSUM)/HSUM 170 D(1,N) = W1*DEL1 + W2*DEL2 171 IF ( PCHST(D(1,N),DEL2) .LE. ZERO) THEN 172 D(1,N) = ZERO 173 ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO) THEN 174C NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. 175 DMAX = THREE*DEL2 176 IF (ABS(D(1,N)) .GT. ABS(DMAX)) D(1,N) = DMAX 177 ENDIF 178C 179C NORMAL RETURN. 180C 181 5000 CONTINUE 182 RETURN 183C------------- LAST LINE OF PCHCI FOLLOWS ------------------------------ 184 END 185