2
SUBROUTINE DPCHCS (SWITCH, N, H, SLOPE, D, INCFD, IERR)
3
C***BEGIN PROLOGUE DPCHCS
5
C***PURPOSE Adjusts derivative values for DPCHIC
6
C***LIBRARY SLATEC (PCHIP)
7
C***TYPE DOUBLE PRECISION (PCHCS-S, DPCHCS-D)
8
C***AUTHOR Fritsch, F. N., (LLNL)
11
C DPCHCS: DPCHIC Monotonicity Switch Derivative Setter.
13
C Called by DPCHIC to adjust the values of D in the vicinity of a
14
C switch in direction of monotonicity, to produce a more "visually
15
C pleasing" curve than that given by DPCHIM .
17
C ----------------------------------------------------------------------
21
C PARAMETER (INCFD = ...)
23
C DOUBLE PRECISION SWITCH, H(N), SLOPE(N), D(INCFD,N)
25
C CALL DPCHCS (SWITCH, N, H, SLOPE, D, INCFD, IERR)
29
C SWITCH -- (input) indicates the amount of control desired over
30
C local excursions from data.
32
C N -- (input) number of data points. (assumes N.GT.2 .)
34
C H -- (input) real*8 array of interval lengths.
35
C SLOPE -- (input) real*8 array of data slopes.
36
C If the data are (X(I),Y(I)), I=1(1)N, then these inputs are:
38
C SLOPE(I) = (Y(I+1)-Y(I))/H(I), I=1(1)N-1.
40
C D -- (input) real*8 array of derivative values at the data points,
41
C as determined by DPCHCI.
42
C (output) derivatives in the vicinity of switches in direction
43
C of monotonicity may be adjusted to produce a more "visually
45
C The value corresponding to X(I) is stored in
46
C D(1+(I-1)*INCFD), I=1(1)N.
47
C No other entries in D are changed.
49
C INCFD -- (input) increment between successive values in D.
50
C This argument is provided primarily for 2-D applications.
52
C IERR -- (output) error flag. should be zero.
53
C If negative, trouble in DPCHSW. (should never happen.)
56
C WARNING: This routine does no validity-checking of arguments.
59
C Fortran intrinsics used: ABS, MAX, MIN.
62
C***ROUTINES CALLED DPCHST, DPCHSW
63
C***REVISION HISTORY (YYMMDD)
65
C 820617 Redesigned to (1) fix problem with lack of continuity
66
C approaching a flat-topped peak (2) be cleaner and
68
C Eliminated subroutines PCHSA and PCHSX in the process.
69
C 820622 1. Limited fact to not exceed one, so computed D is a
70
C convex combination of DPCHCI value and DPCHSD value.
71
C 2. Changed fudge from 1 to 4 (based on experiments).
72
C 820623 Moved PCHSD to an inline function (eliminating MSWTYP).
73
C 820805 Converted to SLATEC library version.
74
C 870707 Corrected conversion to double precision.
75
C 870813 Minor cosmetic changes.
76
C 890411 Added SAVE statements (Vers. 3.2).
77
C 890531 Changed all specific intrinsics to generic. (WRB)
78
C 890831 Modified array declarations. (WRB)
79
C 891006 Modified spacing in computation of DFLOC. (WRB)
80
C 891006 REVISION DATE from Version 3.2
81
C 891214 Prologue converted to Version 4.0 format. (BAB)
82
C 900328 Added TYPE section. (WRB)
83
C 910408 Updated AUTHOR section in prologue. (WRB)
84
C 930503 Improved purpose. (FNF)
85
C***END PROLOGUE DPCHCS
88
C 1. The function DPCHST(ARG1,ARG2) is assumed to return zero if
89
C either argument is zero, +1 if they are of the same sign, and
90
C -1 if they are of opposite sign.
95
INTEGER N, INCFD, IERR
96
DOUBLE PRECISION SWITCH, H(*), SLOPE(*), D(INCFD,*)
98
C DECLARE LOCAL VARIABLES.
100
INTEGER I, INDX, K, NLESS1
101
DOUBLE PRECISION DEL(3), DEXT, DFLOC, DFMX, FACT, FUDGE, ONE,
102
* SLMAX, WTAVE(2), ZERO
103
SAVE ZERO, ONE, FUDGE
104
DOUBLE PRECISION DPCHST
106
C DEFINE INLINE FUNCTION FOR WEIGHTED AVERAGE OF SLOPES.
108
DOUBLE PRECISION DPCHSD, S1, S2, H1, H2
109
DPCHSD(S1,S2,H1,H2) = (H2/(H1+H2))*S1 + (H1/(H1+H2))*S2
113
DATA ZERO /0.D0/, ONE/1.D0/
115
C***FIRST EXECUTABLE STATEMENT DPCHCS
119
C LOOP OVER SEGMENTS.
122
IF ( DPCHST(SLOPE(I-1),SLOPE(I)) ) 100, 300, 900
123
C --------------------------
127
C....... SLOPE SWITCHES MONOTONICITY AT I-TH POINT .....................
129
C DO NOT CHANGE D IF 'UP-DOWN-UP'.
131
IF ( DPCHST(SLOPE(I-2),SLOPE(I)) .GT. ZERO) GO TO 900
132
C --------------------------
134
IF (I .LT. NLESS1) THEN
135
IF ( DPCHST(SLOPE(I+1),SLOPE(I-1)) .GT. ZERO) GO TO 900
136
C ----------------------------
139
C ....... COMPUTE PROVISIONAL VALUE FOR D(1,I).
141
DEXT = DPCHSD (SLOPE(I-1), SLOPE(I), H(I-1), H(I))
143
C ....... DETERMINE WHICH INTERVAL CONTAINS THE EXTREMUM.
145
IF ( DPCHST(DEXT, SLOPE(I-1)) ) 200, 900, 250
146
C -----------------------
149
C DEXT AND SLOPE(I-1) HAVE OPPOSITE SIGNS --
150
C EXTREMUM IS IN (X(I-1),X(I)).
152
C SET UP TO COMPUTE NEW VALUES FOR D(1,I-1) AND D(1,I).
155
* WTAVE(1) = DPCHSD (SLOPE(K-1), SLOPE(K), H(K-1), H(K))
159
C DEXT AND SLOPE(I) HAVE OPPOSITE SIGNS --
160
C EXTREMUM IS IN (X(I),X(I+1)).
162
C SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1).
165
* WTAVE(2) = DPCHSD (SLOPE(K), SLOPE(K+1), H(K), H(K+1))
170
C....... AT LEAST ONE OF SLOPE(I-1) AND SLOPE(I) IS ZERO --
171
C CHECK FOR FLAT-TOPPED PEAK .......................
173
IF (I .EQ. NLESS1) GO TO 900
174
IF ( DPCHST(SLOPE(I-1), SLOPE(I+1)) .GE. ZERO) GO TO 900
175
C -----------------------------
177
C WE HAVE FLAT-TOPPED PEAK ON (X(I),X(I+1)).
179
C SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1).
180
WTAVE(1) = DPCHSD (SLOPE(K-1), SLOPE(K), H(K-1), H(K))
181
WTAVE(2) = DPCHSD (SLOPE(K), SLOPE(K+1), H(K), H(K+1))
185
C....... AT THIS POINT WE HAVE DETERMINED THAT THERE WILL BE AN EXTREMUM
186
C ON (X(K),X(K+1)), WHERE K=I OR I-1, AND HAVE SET ARRAY WTAVE--
187
C WTAVE(1) IS A WEIGHTED AVERAGE OF SLOPE(K-1) AND SLOPE(K),
189
C WTAVE(2) IS A WEIGHTED AVERAGE OF SLOPE(K) AND SLOPE(K+1),
192
SLMAX = ABS(SLOPE(K))
193
IF (K .GT. 1) SLMAX = MAX( SLMAX, ABS(SLOPE(K-1)) )
194
IF (K.LT.NLESS1) SLMAX = MAX( SLMAX, ABS(SLOPE(K+1)) )
196
IF (K .GT. 1) DEL(1) = SLOPE(K-1) / SLMAX
197
DEL(2) = SLOPE(K) / SLMAX
198
IF (K.LT.NLESS1) DEL(3) = SLOPE(K+1) / SLMAX
200
IF ((K.GT.1) .AND. (K.LT.NLESS1)) THEN
201
C NORMAL CASE -- EXTREMUM IS NOT IN A BOUNDARY INTERVAL.
202
FACT = FUDGE* ABS(DEL(3)*(DEL(1)-DEL(2))*(WTAVE(2)/SLMAX))
203
D(1,K) = D(1,K) + MIN(FACT,ONE)*(WTAVE(1) - D(1,K))
204
FACT = FUDGE* ABS(DEL(1)*(DEL(3)-DEL(2))*(WTAVE(1)/SLMAX))
205
D(1,K+1) = D(1,K+1) + MIN(FACT,ONE)*(WTAVE(2) - D(1,K+1))
207
C SPECIAL CASE K=1 (WHICH CAN OCCUR ONLY IF I=2) OR
208
C K=NLESS1 (WHICH CAN OCCUR ONLY IF I=NLESS1).
209
FACT = FUDGE* ABS(DEL(2))
210
D(1,I) = MIN(FACT,ONE) * WTAVE(I-K+1)
211
C NOTE THAT I-K+1 = 1 IF K=I (=NLESS1),
212
C I-K+1 = 2 IF K=I-1(=1).
216
C....... ADJUST IF NECESSARY TO LIMIT EXCURSIONS FROM DATA.
218
IF (SWITCH .LE. ZERO) GO TO 900
220
DFLOC = H(K)*ABS(SLOPE(K))
221
IF (K .GT. 1) DFLOC = MAX( DFLOC, H(K-1)*ABS(SLOPE(K-1)) )
222
IF (K.LT.NLESS1) DFLOC = MAX( DFLOC, H(K+1)*ABS(SLOPE(K+1)) )
225
C INDX = 1 IF K=I, 2 IF K=I-1.
226
C ---------------------------------------------------------------
227
CALL DPCHSW(DFMX, INDX, D(1,K), D(1,K+1), H(K), SLOPE(K), IERR)
228
C ---------------------------------------------------------------
229
IF (IERR .NE. 0) RETURN
231
C....... END OF SEGMENT LOOP.
236
C------------- LAST LINE OF DPCHCS FOLLOWS -----------------------------
added
removed
Lib/Slatec/slatec/dpchcs.f
