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
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
|
SUBROUTINE dinvr(status,x,fx,qleft,qhi)
C**********************************************************************
C
C SUBROUTINE DINVR(STATUS, X, FX, QLEFT, QHI)
C Double precision
C bounds the zero of the function and invokes zror
C Reverse Communication
C
C
C Function
C
C
C Bounds the function and invokes ZROR to perform the zero
C finding. STINVR must have been called before this routine
C in order to set its parameters.
C
C
C Arguments
C
C
C STATUS <--> At the beginning of a zero finding problem, STATUS
C should be set to 0 and INVR invoked. (The value
C of parameters other than X will be ignored on this cal
C
C When INVR needs the function evaluated, it will set
C STATUS to 1 and return. The value of the function
C should be set in FX and INVR again called without
C changing any of its other parameters.
C
C When INVR has finished without error, it will return
C with STATUS 0. In that case X is approximately a root
C of F(X).
C
C If INVR cannot bound the function, it returns status
C -1 and sets QLEFT and QHI.
C INTEGER STATUS
C
C X <-- The value of X at which F(X) is to be evaluated.
C DOUBLE PRECISION X
C
C FX --> The value of F(X) calculated when INVR returns with
C STATUS = 1.
C DOUBLE PRECISION FX
C
C QLEFT <-- Defined only if QMFINV returns .FALSE. In that
C case it is .TRUE. If the stepping search terminated
C unsucessfully at SMALL. If it is .FALSE. the search
C terminated unsucessfully at BIG.
C QLEFT is LOGICAL
C
C QHI <-- Defined only if QMFINV returns .FALSE. In that
C case it is .TRUE. if F(X) .GT. Y at the termination
C of the search and .FALSE. if F(X) .LT. Y at the
C termination of the search.
C QHI is LOGICAL
C
C**********************************************************************
C Modified by S. Steer INRIA 1998,to replace ASSIGN instruction by
c Computed GOTO
C**********************************************************************
include '../stack.h'
C .. Scalar Arguments ..
DOUBLE PRECISION fx,x,zabsst,zabsto,zbig,zrelst,zrelto,zsmall,
+ zstpmu
INTEGER status
LOGICAL qhi,qleft
C ..
C .. Local Scalars ..
DOUBLE PRECISION absstp,abstol,big,fbig,fsmall,relstp,reltol,
+ small,step,stpmul,xhi,xlb,xlo,xsave,xub,yy,zx,zy,
+ zz
INTEGER i99999
LOGICAL qbdd,qcond,qdum1,qdum2,qincr,qlim,qok,qup
C ..
C .. External Subroutines ..
EXTERNAL dstzr,dzror
C ..
C .. Intrinsic Functions ..
INTRINSIC abs,max,min
C ..
C .. Statement Functions ..
LOGICAL qxmon
C ..
C .. Save statement ..
SAVE
C ..
C .. Statement Function definitions ..
qxmon(zx,zy,zz) = zx .LE. zy .AND. zy .LE. zz
C ..
C .. Executable Statements ..
IF (status.GT.0) GO TO 310
qcond = .NOT. qxmon(small,x,big)
IF (qcond) then
call basout(io,wte,' SMALL, X, BIG not monotone in INVR')
status = -100
return
endif
xsave = x
C
C See that SMALL and BIG bound the zero and set QINCR
C
x = small
C GET-FUNCTION-VALUE
c ASSIGN 10 TO i99999
i99999=1
GO TO 300
10 fsmall = fx
x = big
C GET-FUNCTION-VALUE
c ASSIGN 20 TO i99999
i99999=2
GO TO 300
20 fbig = fx
qincr = fbig .GT. fsmall
IF (.NOT. (qincr)) GO TO 50
IF (fsmall.LE.0.0D0) GO TO 30
status = -1
qleft = .TRUE.
qhi = .TRUE.
RETURN
30 IF (fbig.GE.0.0D0) GO TO 40
status = -1
qleft = .FALSE.
qhi = .FALSE.
RETURN
40 GO TO 80
50 IF (fsmall.GE.0.0D0) GO TO 60
status = -1
qleft = .TRUE.
qhi = .FALSE.
RETURN
60 IF (fbig.LE.0.0D0) GO TO 70
status = -1
qleft = .FALSE.
qhi = .TRUE.
RETURN
70 CONTINUE
80 x = xsave
step = max(absstp,relstp*abs(x))
C YY = F(X) - Y
C GET-FUNCTION-VALUE
c ASSIGN 90 TO i99999
i99999=3
GO TO 300
90 yy = fx
IF (.NOT. (yy.EQ.0.0D0)) GO TO 100
status = 0
qok = .TRUE.
RETURN
100 qup = (qincr .AND. (yy.LT.0.0D0)) .OR.
+ (.NOT.qincr .AND. (yy.GT.0.0D0))
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C HANDLE CASE IN WHICH WE MUST STEP HIGHER
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
IF (.NOT. (qup)) GO TO 170
xlb = xsave
xub = min(xlb+step,big)
GO TO 120
110 IF (qcond) GO TO 150
C YY = F(XUB) - Y
120 x = xub
C GET-FUNCTION-VALUE
c ASSIGN 130 TO i99999
i99999=4
GO TO 300
130 yy = fx
qbdd = (qincr .AND. (yy.GE.0.0D0)) .OR.
+ (.NOT.qincr .AND. (yy.LE.0.0D0))
qlim = xub .GE. big
qcond = qbdd .OR. qlim
IF (qcond) GO TO 140
step = stpmul*step
xlb = xub
xub = min(xlb+step,big)
140 GO TO 110
150 IF (.NOT. (qlim.AND..NOT.qbdd)) GO TO 160
status = -1
qleft = .FALSE.
qhi = .NOT. qincr
x = big
RETURN
160 GO TO 240
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C HANDLE CASE IN WHICH WE MUST STEP LOWER
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
170 xub = xsave
xlb = max(xub-step,small)
GO TO 190
180 IF (qcond) GO TO 220
C YY = F(XLB) - Y
190 x = xlb
C GET-FUNCTION-VALUE
c ASSIGN 200 TO i99999
i99999=5
GO TO 300
200 yy = fx
qbdd = (qincr .AND. (yy.LE.0.0D0)) .OR.
+ (.NOT.qincr .AND. (yy.GE.0.0D0))
qlim = xlb .LE. small
qcond = qbdd .OR. qlim
IF (qcond) GO TO 210
step = stpmul*step
xub = xlb
xlb = max(xub-step,small)
210 GO TO 180
220 IF (.NOT. (qlim.AND..NOT.qbdd)) GO TO 230
status = -1
qleft = .TRUE.
qhi = qincr
x = small
RETURN
230 CONTINUE
240 CALL dstzr(xlb,xub,abstol,reltol)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C IF WE REACH HERE, XLB AND XUB BOUND THE ZERO OF F.
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
status = 0
GO TO 260
250 IF (.NOT. (status.EQ.1)) GO TO 290
260 CALL dzror(status,x,fx,xlo,xhi,qdum1,qdum2)
IF (.NOT. (status.EQ.1)) GO TO 280
C GET-FUNCTION-VALUE
c ASSIGN 270 TO i99999
i99999=6
GO TO 300
270 CONTINUE
280 GO TO 250
290 x = xlo
status = 0
RETURN
ENTRY dstinv(zsmall,zbig,zabsst,zrelst,zstpmu,zabsto,zrelto)
C**********************************************************************
C
C SUBROUTINE DSTINV( SMALL, BIG, ABSSTP, RELSTP, STPMUL,
C + ABSTOL, RELTOL )
C Double Precision - SeT INverse finder - Reverse Communication
C
C
C Function
C
C
C Concise Description - Given a monotone function F finds X
C such that F(X) = Y. Uses Reverse communication -- see invr.
C This routine sets quantities needed by INVR.
C
C More Precise Description of INVR -
C
C F must be a monotone function, the results of QMFINV are
C otherwise undefined. QINCR must be .TRUE. if F is non-
C decreasing and .FALSE. if F is non-increasing.
C
C QMFINV will return .TRUE. if and only if F(SMALL) and
C F(BIG) bracket Y, i. e.,
C QINCR is .TRUE. and F(SMALL).LE.Y.LE.F(BIG) or
C QINCR is .FALSE. and F(BIG).LE.Y.LE.F(SMALL)
C
C if QMFINV returns .TRUE., then the X returned satisfies
C the following condition. let
C TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
C then if QINCR is .TRUE.,
C F(X-TOL(X)) .LE. Y .LE. F(X+TOL(X))
C and if QINCR is .FALSE.
C F(X-TOL(X)) .GE. Y .GE. F(X+TOL(X))
C
C
C Arguments
C
C
C SMALL --> The left endpoint of the interval to be
C searched for a solution.
C SMALL is DOUBLE PRECISION
C
C BIG --> The right endpoint of the interval to be
C searched for a solution.
C BIG is DOUBLE PRECISION
C
C ABSSTP, RELSTP --> The initial step size in the search
C is MAX(ABSSTP,RELSTP*ABS(X)). See algorithm.
C ABSSTP is DOUBLE PRECISION
C RELSTP is DOUBLE PRECISION
C
C STPMUL --> When a step doesn't bound the zero, the step
C size is multiplied by STPMUL and another step
C taken. A popular value is 2.0
C DOUBLE PRECISION STPMUL
C
C ABSTOL, RELTOL --> Two numbers that determine the accuracy
C of the solution. See function for a precise definition.
C ABSTOL is DOUBLE PRECISION
C RELTOL is DOUBLE PRECISION
C
C
C Method
C
C
C Compares F(X) with Y for the input value of X then uses QINCR
C to determine whether to step left or right to bound the
C desired x. the initial step size is
C MAX(ABSSTP,RELSTP*ABS(S)) for the input value of X.
C Iteratively steps right or left until it bounds X.
C At each step which doesn't bound X, the step size is doubled.
C The routine is careful never to step beyond SMALL or BIG. If
C it hasn't bounded X at SMALL or BIG, QMFINV returns .FALSE.
C after setting QLEFT and QHI.
C
C If X is successfully bounded then Algorithm R of the paper
C 'Two Efficient Algorithms with Guaranteed Convergence for
C Finding a Zero of a Function' by J. C. P. Bus and
C T. J. Dekker in ACM Transactions on Mathematical
C Software, Volume 1, No. 4 page 330 (DEC. '75) is employed
C to find the zero of the function F(X)-Y. This is routine
C QRZERO.
C
C**********************************************************************
small = zsmall
big = zbig
absstp = zabsst
relstp = zrelst
stpmul = zstpmu
abstol = zabsto
reltol = zrelto
RETURN
C(jpc) STOP '*** EXECUTION FLOWING INTO FLECS PROCEDURES ***'
C TO GET-FUNCTION-VALUE
300 status = 1
RETURN
310 CONTINUE
goto(10,20,90,130,200,270) i99999
c GO TO i99999
END
|
