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