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
|
*DECK PCHSW
SUBROUTINE PCHSW (DFMAX, IEXTRM, D1, D2, H, SLOPE, IERR)
C***BEGIN PROLOGUE PCHSW
C***SUBSIDIARY
C***PURPOSE Limits excursion from data for PCHCS
C***LIBRARY SLATEC (PCHIP)
C***TYPE SINGLE PRECISION (PCHSW-S, DPCHSW-D)
C***AUTHOR Fritsch, F. N., (LLNL)
C***DESCRIPTION
C
C PCHSW: PCHCS Switch Excursion Limiter.
C
C Called by PCHCS to adjust D1 and D2 if necessary to insure that
C the extremum on this interval is not further than DFMAX from the
C extreme data value.
C
C ----------------------------------------------------------------------
C
C Calling sequence:
C
C INTEGER IEXTRM, IERR
C REAL DFMAX, D1, D2, H, SLOPE
C
C CALL PCHSW (DFMAX, IEXTRM, D1, D2, H, SLOPE, IERR)
C
C Parameters:
C
C DFMAX -- (input) maximum allowed difference between F(IEXTRM) and
C the cubic determined by derivative values D1,D2. (assumes
C DFMAX.GT.0.)
C
C IEXTRM -- (input) index of the extreme data value. (assumes
C IEXTRM = 1 or 2 . Any value .NE.1 is treated as 2.)
C
C D1,D2 -- (input) derivative values at the ends of the interval.
C (Assumes D1*D2 .LE. 0.)
C (output) may be modified if necessary to meet the restriction
C imposed by DFMAX.
C
C H -- (input) interval length. (Assumes H.GT.0.)
C
C SLOPE -- (input) data slope on the interval.
C
C IERR -- (output) error flag. should be zero.
C If IERR=-1, assumption on D1 and D2 is not satisfied.
C If IERR=-2, quadratic equation locating extremum has
C negative discriminant (should never occur).
C
C -------
C WARNING: This routine does no validity-checking of arguments.
C -------
C
C Fortran intrinsics used: ABS, SIGN, SQRT.
C
C***SEE ALSO PCHCS
C***ROUTINES CALLED R1MACH, XERMSG
C***REVISION HISTORY (YYMMDD)
C 820218 DATE WRITTEN
C 820805 Converted to SLATEC library version.
C 870707 Replaced DATA statement for SMALL with a use of R1MACH.
C 890411 1. Added SAVE statements (Vers. 3.2).
C 2. Added REAL R1MACH for consistency with D.P. version.
C 890411 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
C 900328 Added TYPE section. (WRB)
C 910408 Updated AUTHOR and DATE WRITTEN sections in prologue. (WRB)
C 920526 Eliminated possible divide by zero problem. (FNF)
C 930503 Improved purpose. (FNF)
C***END PROLOGUE PCHSW
C
C**End
C
C DECLARE ARGUMENTS.
C
INTEGER IEXTRM, IERR
REAL DFMAX, D1, D2, H, SLOPE
C
C DECLARE LOCAL VARIABLES.
C
REAL CP, FACT, HPHI, LAMBDA, NU, ONE, PHI, RADCAL, RHO, SIGMA,
* SMALL, THAT, THIRD, THREE, TWO, ZERO
SAVE ZERO, ONE, TWO, THREE, FACT
SAVE THIRD
REAL R1MACH
C
DATA ZERO /0./, ONE /1./, TWO /2./, THREE /3./, FACT /100./
C THIRD SHOULD BE SLIGHTLY LESS THAN 1/3.
DATA THIRD /0.33333/
C
C NOTATION AND GENERAL REMARKS.
C
C RHO IS THE RATIO OF THE DATA SLOPE TO THE DERIVATIVE BEING TESTED.
C LAMBDA IS THE RATIO OF D2 TO D1.
C THAT = T-HAT(RHO) IS THE NORMALIZED LOCATION OF THE EXTREMUM.
C PHI IS THE NORMALIZED VALUE OF P(X)-F1 AT X = XHAT = X-HAT(RHO),
C WHERE THAT = (XHAT - X1)/H .
C THAT IS, P(XHAT)-F1 = D*H*PHI, WHERE D=D1 OR D2.
C SIMILARLY, P(XHAT)-F2 = D*H*(PHI-RHO) .
C
C SMALL SHOULD BE A FEW ORDERS OF MAGNITUDE GREATER THAN MACHEPS.
C***FIRST EXECUTABLE STATEMENT PCHSW
SMALL = FACT*R1MACH(4)
C
C DO MAIN CALCULATION.
C
IF (D1 .EQ. ZERO) THEN
C
C SPECIAL CASE -- D1.EQ.ZERO .
C
C IF D2 IS ALSO ZERO, THIS ROUTINE SHOULD NOT HAVE BEEN CALLED.
IF (D2 .EQ. ZERO) GO TO 5001
C
RHO = SLOPE/D2
C EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 .
IF (RHO .GE. THIRD) GO TO 5000
THAT = (TWO*(THREE*RHO-ONE)) / (THREE*(TWO*RHO-ONE))
PHI = THAT**2 * ((THREE*RHO-ONE)/THREE)
C
C CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 .
IF (IEXTRM .NE. 1) PHI = PHI - RHO
C
C TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY.
HPHI = H * ABS(PHI)
IF (HPHI*ABS(D2) .GT. DFMAX) THEN
C AT THIS POINT, HPHI.GT.0, SO DIVIDE IS OK.
D2 = SIGN (DFMAX/HPHI, D2)
ENDIF
ELSE
C
RHO = SLOPE/D1
LAMBDA = -D2/D1
IF (D2 .EQ. ZERO) THEN
C
C SPECIAL CASE -- D2.EQ.ZERO .
C
C EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 .
IF (RHO .GE. THIRD) GO TO 5000
CP = TWO - THREE*RHO
NU = ONE - TWO*RHO
THAT = ONE / (THREE*NU)
ELSE
IF (LAMBDA .LE. ZERO) GO TO 5001
C
C NORMAL CASE -- D1 AND D2 BOTH NONZERO, OPPOSITE SIGNS.
C
NU = ONE - LAMBDA - TWO*RHO
SIGMA = ONE - RHO
CP = NU + SIGMA
IF (ABS(NU) .GT. SMALL) THEN
RADCAL = (NU - (TWO*RHO+ONE))*NU + SIGMA**2
IF (RADCAL .LT. ZERO) GO TO 5002
THAT = (CP - SQRT(RADCAL)) / (THREE*NU)
ELSE
THAT = ONE/(TWO*SIGMA)
ENDIF
ENDIF
PHI = THAT*((NU*THAT - CP)*THAT + ONE)
C
C CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 .
IF (IEXTRM .NE. 1) PHI = PHI - RHO
C
C TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY.
HPHI = H * ABS(PHI)
IF (HPHI*ABS(D1) .GT. DFMAX) THEN
C AT THIS POINT, HPHI.GT.0, SO DIVIDE IS OK.
D1 = SIGN (DFMAX/HPHI, D1)
D2 = -LAMBDA*D1
ENDIF
ENDIF
C
C NORMAL RETURN.
C
5000 CONTINUE
IERR = 0
RETURN
C
C ERROR RETURNS.
C
5001 CONTINUE
C D1 AND D2 BOTH ZERO, OR BOTH NONZERO AND SAME SIGN.
IERR = -1
CALL XERMSG ('SLATEC', 'PCHSW', 'D1 AND/OR D2 INVALID', IERR, 1)
RETURN
C
5002 CONTINUE
C NEGATIVE VALUE OF RADICAL (SHOULD NEVER OCCUR).
IERR = -2
CALL XERMSG ('SLATEC', 'PCHSW', 'NEGATIVE RADICAL', IERR, 1)
RETURN
C------------- LAST LINE OF PCHSW FOLLOWS ------------------------------
END
|
