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
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
|
*DECK PCHCS
SUBROUTINE PCHCS (SWITCH, N, H, SLOPE, D, INCFD, IERR)
C***BEGIN PROLOGUE PCHCS
C***SUBSIDIARY
C***PURPOSE Adjusts derivative values for PCHIC
C***LIBRARY SLATEC (PCHIP)
C***TYPE SINGLE PRECISION (PCHCS-S, DPCHCS-D)
C***AUTHOR Fritsch, F. N., (LLNL)
C***DESCRIPTION
C
C PCHCS: PCHIC Monotonicity Switch Derivative Setter.
C
C Called by PCHIC to adjust the values of D in the vicinity of a
C switch in direction of monotonicity, to produce a more "visually
C pleasing" curve than that given by PCHIM .
C
C ----------------------------------------------------------------------
C
C Calling sequence:
C
C PARAMETER (INCFD = ...)
C INTEGER N, IERR
C REAL SWITCH, H(N), SLOPE(N), D(INCFD,N)
C
C CALL PCHCS (SWITCH, N, H, SLOPE, D, INCFD, IERR)
C
C Parameters:
C
C SWITCH -- (input) indicates the amount of control desired over
C local excursions from data.
C
C N -- (input) number of data points. (assumes N.GT.2 .)
C
C H -- (input) real array of interval lengths.
C SLOPE -- (input) real 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 -- (input) real array of derivative values at the data points,
C as determined by PCHCI.
C (output) derivatives in the vicinity of switches in direction
C of monotonicity may be adjusted to produce a more "visually
C pleasing" curve.
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 IERR -- (output) error flag. should be zero.
C If negative, trouble in PCHSW. (should never happen.)
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 PCHIC
C***ROUTINES CALLED PCHST, PCHSW
C***REVISION HISTORY (YYMMDD)
C 820218 DATE WRITTEN
C 820617 Redesigned to (1) fix problem with lack of continuity
C approaching a flat-topped peak (2) be cleaner and
C easier to verify.
C Eliminated subroutines PCHSA and PCHSX in the process.
C 820622 1. Limited fact to not exceed one, so computed D is a
C convex combination of PCHCI value and PCHSD value.
C 2. Changed fudge from 1 to 4 (based on experiments).
C 820623 Moved PCHSD to an inline function (eliminating MSWTYP).
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 PCHCS
C
C Programming notes:
C 1. The function PCHST(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, IERR
REAL SWITCH, H(*), SLOPE(*), D(INCFD,*)
C
C DECLARE LOCAL VARIABLES.
C
INTEGER I, INDX, K, NLESS1
REAL DEL(3), DEXT, DFLOC, DFMX, FACT, FUDGE, ONE, SLMAX,
* WTAVE(2), ZERO
SAVE ZERO, ONE, FUDGE
REAL PCHST
C
C DEFINE INLINE FUNCTION FOR WEIGHTED AVERAGE OF SLOPES.
C
REAL PCHSD, S1, S2, H1, H2
PCHSD(S1,S2,H1,H2) = (H2/(H1+H2))*S1 + (H1/(H1+H2))*S2
C
C INITIALIZE.
C
DATA ZERO /0./, ONE /1./
DATA FUDGE /4./
C***FIRST EXECUTABLE STATEMENT PCHCS
IERR = 0
NLESS1 = N - 1
C
C LOOP OVER SEGMENTS.
C
DO 900 I = 2, NLESS1
IF ( PCHST(SLOPE(I-1),SLOPE(I)) ) 100, 300, 900
C --------------------------
C
100 CONTINUE
C
C....... SLOPE SWITCHES MONOTONICITY AT I-TH POINT .....................
C
C DO NOT CHANGE D IF 'UP-DOWN-UP'.
IF (I .GT. 2) THEN
IF ( PCHST(SLOPE(I-2),SLOPE(I)) .GT. ZERO) GO TO 900
C --------------------------
ENDIF
IF (I .LT. NLESS1) THEN
IF ( PCHST(SLOPE(I+1),SLOPE(I-1)) .GT. ZERO) GO TO 900
C ----------------------------
ENDIF
C
C ....... COMPUTE PROVISIONAL VALUE FOR D(1,I).
C
DEXT = PCHSD (SLOPE(I-1), SLOPE(I), H(I-1), H(I))
C
C ....... DETERMINE WHICH INTERVAL CONTAINS THE EXTREMUM.
C
IF ( PCHST(DEXT, SLOPE(I-1)) ) 200, 900, 250
C -----------------------
C
200 CONTINUE
C DEXT AND SLOPE(I-1) HAVE OPPOSITE SIGNS --
C EXTREMUM IS IN (X(I-1),X(I)).
K = I-1
C SET UP TO COMPUTE NEW VALUES FOR D(1,I-1) AND D(1,I).
WTAVE(2) = DEXT
IF (K .GT. 1)
* WTAVE(1) = PCHSD (SLOPE(K-1), SLOPE(K), H(K-1), H(K))
GO TO 400
C
250 CONTINUE
C DEXT AND SLOPE(I) HAVE OPPOSITE SIGNS --
C EXTREMUM IS IN (X(I),X(I+1)).
K = I
C SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1).
WTAVE(1) = DEXT
IF (K .LT. NLESS1)
* WTAVE(2) = PCHSD (SLOPE(K), SLOPE(K+1), H(K), H(K+1))
GO TO 400
C
300 CONTINUE
C
C....... AT LEAST ONE OF SLOPE(I-1) AND SLOPE(I) IS ZERO --
C CHECK FOR FLAT-TOPPED PEAK .......................
C
IF (I .EQ. NLESS1) GO TO 900
IF ( PCHST(SLOPE(I-1), SLOPE(I+1)) .GE. ZERO) GO TO 900
C -----------------------------
C
C WE HAVE FLAT-TOPPED PEAK ON (X(I),X(I+1)).
K = I
C SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1).
WTAVE(1) = PCHSD (SLOPE(K-1), SLOPE(K), H(K-1), H(K))
WTAVE(2) = PCHSD (SLOPE(K), SLOPE(K+1), H(K), H(K+1))
C
400 CONTINUE
C
C....... AT THIS POINT WE HAVE DETERMINED THAT THERE WILL BE AN EXTREMUM
C ON (X(K),X(K+1)), WHERE K=I OR I-1, AND HAVE SET ARRAY WTAVE--
C WTAVE(1) IS A WEIGHTED AVERAGE OF SLOPE(K-1) AND SLOPE(K),
C IF K.GT.1
C WTAVE(2) IS A WEIGHTED AVERAGE OF SLOPE(K) AND SLOPE(K+1),
C IF K.LT.N-1
C
SLMAX = ABS(SLOPE(K))
IF (K .GT. 1) SLMAX = MAX( SLMAX, ABS(SLOPE(K-1)) )
IF (K.LT.NLESS1) SLMAX = MAX( SLMAX, ABS(SLOPE(K+1)) )
C
IF (K .GT. 1) DEL(1) = SLOPE(K-1) / SLMAX
DEL(2) = SLOPE(K) / SLMAX
IF (K.LT.NLESS1) DEL(3) = SLOPE(K+1) / SLMAX
C
IF ((K.GT.1) .AND. (K.LT.NLESS1)) THEN
C NORMAL CASE -- EXTREMUM IS NOT IN A BOUNDARY INTERVAL.
FACT = FUDGE* ABS(DEL(3)*(DEL(1)-DEL(2))*(WTAVE(2)/SLMAX))
D(1,K) = D(1,K) + MIN(FACT,ONE)*(WTAVE(1) - D(1,K))
FACT = FUDGE* ABS(DEL(1)*(DEL(3)-DEL(2))*(WTAVE(1)/SLMAX))
D(1,K+1) = D(1,K+1) + MIN(FACT,ONE)*(WTAVE(2) - D(1,K+1))
ELSE
C SPECIAL CASE K=1 (WHICH CAN OCCUR ONLY IF I=2) OR
C K=NLESS1 (WHICH CAN OCCUR ONLY IF I=NLESS1).
FACT = FUDGE* ABS(DEL(2))
D(1,I) = MIN(FACT,ONE) * WTAVE(I-K+1)
C NOTE THAT I-K+1 = 1 IF K=I (=NLESS1),
C I-K+1 = 2 IF K=I-1(=1).
ENDIF
C
C
C....... ADJUST IF NECESSARY TO LIMIT EXCURSIONS FROM DATA.
C
IF (SWITCH .LE. ZERO) GO TO 900
C
DFLOC = H(K)*ABS(SLOPE(K))
IF (K .GT. 1) DFLOC = MAX( DFLOC, H(K-1)*ABS(SLOPE(K-1)) )
IF (K.LT.NLESS1) DFLOC = MAX( DFLOC, H(K+1)*ABS(SLOPE(K+1)) )
DFMX = SWITCH*DFLOC
INDX = I-K+1
C INDX = 1 IF K=I, 2 IF K=I-1.
C ---------------------------------------------------------------
CALL PCHSW (DFMX, INDX, D(1,K), D(1,K+1), H(K), SLOPE(K), IERR)
C ---------------------------------------------------------------
IF (IERR .NE. 0) RETURN
C
C....... END OF SEGMENT LOOP.
C
900 CONTINUE
C
RETURN
C------------- LAST LINE OF PCHCS FOLLOWS ------------------------------
END
|
