1
SUBROUTINE dinvr(status,x,fx,qleft,qhi)
2
C**********************************************************************
4
C SUBROUTINE DINVR(STATUS, X, FX, QLEFT, QHI)
6
C bounds the zero of the function and invokes zror
7
C Reverse Communication
13
C Bounds the function and invokes ZROR to perform the zero
14
C finding. STINVR must have been called before this routine
15
C in order to set its parameters.
21
C STATUS <--> At the beginning of a zero finding problem, STATUS
22
C should be set to 0 and INVR invoked. (The value
23
C of parameters other than X will be ignored on this cal
25
C When INVR needs the function evaluated, it will set
26
C STATUS to 1 and return. The value of the function
27
C should be set in FX and INVR again called without
28
C changing any of its other parameters.
30
C When INVR has finished without error, it will return
31
C with STATUS 0. In that case X is approximately a root
34
C If INVR cannot bound the function, it returns status
35
C -1 and sets QLEFT and QHI.
38
C X <-- The value of X at which F(X) is to be evaluated.
41
C FX --> The value of F(X) calculated when INVR returns with
45
C QLEFT <-- Defined only if QMFINV returns .FALSE. In that
46
C case it is .TRUE. If the stepping search terminated
47
C unsucessfully at SMALL. If it is .FALSE. the search
48
C terminated unsucessfully at BIG.
51
C QHI <-- Defined only if QMFINV returns .FALSE. In that
52
C case it is .TRUE. if F(X) .GT. Y at the termination
53
C of the search and .FALSE. if F(X) .LT. Y at the
54
C termination of the search.
58
C**********************************************************************
59
C .. Scalar Arguments ..
60
DOUBLE PRECISION fx,x,zabsst,zabsto,zbig,zrelst,zrelto,zsmall,
66
DOUBLE PRECISION absstp,abstol,big,fbig,fsmall,relstp,reltol,
67
+ small,step,stpmul,xhi,xlb,xlo,xsave,xub,yy,zx,zy,
70
LOGICAL qbdd,qcond,qdum1,qdum2,qincr,qlim,qok,qup
72
C .. External Subroutines ..
75
C .. Intrinsic Functions ..
78
C .. Statement Functions ..
81
C .. Save statement ..
84
C .. Statement Function definitions ..
85
qxmon(zx,zy,zz) = zx .LE. zy .AND. zy .LE. zz
87
C .. Executable Statements ..
89
IF (status.GT.0) GO TO 310
91
qcond = .NOT. qxmon(small,x,big)
92
IF (qcond) STOP ' SMALL, X, BIG not monotone in INVR'
95
C See that SMALL and BIG bound the zero and set QINCR
109
qincr = fbig .GT. fsmall
110
IF (.NOT. (qincr)) GO TO 50
111
IF (fsmall.LE.0.0D0) GO TO 30
117
30 IF (fbig.GE.0.0D0) GO TO 40
125
50 IF (fsmall.GE.0.0D0) GO TO 60
131
60 IF (fbig.LE.0.0D0) GO TO 70
139
step = max(absstp,relstp*abs(x))
146
IF (.NOT. (yy.EQ.0.0D0)) GO TO 100
151
100 qup = (qincr .AND. (yy.LT.0.0D0)) .OR.
152
+ (.NOT.qincr .AND. (yy.GT.0.0D0))
153
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
155
C HANDLE CASE IN WHICH WE MUST STEP HIGHER
157
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
158
IF (.NOT. (qup)) GO TO 170
160
xub = min(xlb+step,big)
163
110 IF (qcond) GO TO 150
171
qbdd = (qincr .AND. (yy.GE.0.0D0)) .OR.
172
+ (.NOT.qincr .AND. (yy.LE.0.0D0))
174
qcond = qbdd .OR. qlim
178
xub = min(xlb+step,big)
181
150 IF (.NOT. (qlim.AND..NOT.qbdd)) GO TO 160
189
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
191
C HANDLE CASE IN WHICH WE MUST STEP LOWER
193
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
195
xlb = max(xub-step,small)
198
180 IF (qcond) GO TO 220
206
qbdd = (qincr .AND. (yy.LE.0.0D0)) .OR.
207
+ (.NOT.qincr .AND. (yy.GE.0.0D0))
208
qlim = xlb .LE. small
209
qcond = qbdd .OR. qlim
213
xlb = max(xub-step,small)
216
220 IF (.NOT. (qlim.AND..NOT.qbdd)) GO TO 230
224
240 CALL dstzr(xlb,xub,abstol,reltol)
225
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
227
C IF WE REACH HERE, XLB AND XUB BOUND THE ZERO OF F.
229
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
233
250 IF (.NOT. (status.EQ.1)) GO TO 290
234
260 CALL dzror(status,x,fx,xlo,xhi,qdum1,qdum2)
235
IF (.NOT. (status.EQ.1)) GO TO 280
247
ENTRY dstinv(zsmall,zbig,zabsst,zrelst,zstpmu,zabsto,zrelto)
248
C**********************************************************************
250
C SUBROUTINE DSTINV( SMALL, BIG, ABSSTP, RELSTP, STPMUL,
252
C Double Precision - SeT INverse finder - Reverse Communication
258
C Concise Description - Given a monotone function F finds X
259
C such that F(X) = Y. Uses Reverse communication -- see invr.
260
C This routine sets quantities needed by INVR.
262
C More Precise Description of INVR -
264
C F must be a monotone function, the results of QMFINV are
265
C otherwise undefined. QINCR must be .TRUE. if F is non-
266
C decreasing and .FALSE. if F is non-increasing.
268
C QMFINV will return .TRUE. if and only if F(SMALL) and
269
C F(BIG) bracket Y, i. e.,
270
C QINCR is .TRUE. and F(SMALL).LE.Y.LE.F(BIG) or
271
C QINCR is .FALSE. and F(BIG).LE.Y.LE.F(SMALL)
273
C if QMFINV returns .TRUE., then the X returned satisfies
274
C the following condition. let
275
C TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
276
C then if QINCR is .TRUE.,
277
C F(X-TOL(X)) .LE. Y .LE. F(X+TOL(X))
278
C and if QINCR is .FALSE.
279
C F(X-TOL(X)) .GE. Y .GE. F(X+TOL(X))
285
C SMALL --> The left endpoint of the interval to be
286
C searched for a solution.
287
C SMALL is DOUBLE PRECISION
289
C BIG --> The right endpoint of the interval to be
290
C searched for a solution.
291
C BIG is DOUBLE PRECISION
293
C ABSSTP, RELSTP --> The initial step size in the search
294
C is MAX(ABSSTP,RELSTP*ABS(X)). See algorithm.
295
C ABSSTP is DOUBLE PRECISION
296
C RELSTP is DOUBLE PRECISION
298
C STPMUL --> When a step doesn't bound the zero, the step
299
C size is multiplied by STPMUL and another step
300
C taken. A popular value is 2.0
301
C DOUBLE PRECISION STPMUL
303
C ABSTOL, RELTOL --> Two numbers that determine the accuracy
304
C of the solution. See function for a precise definition.
305
C ABSTOL is DOUBLE PRECISION
306
C RELTOL is DOUBLE PRECISION
312
C Compares F(X) with Y for the input value of X then uses QINCR
313
C to determine whether to step left or right to bound the
314
C desired x. the initial step size is
315
C MAX(ABSSTP,RELSTP*ABS(S)) for the input value of X.
316
C Iteratively steps right or left until it bounds X.
317
C At each step which doesn't bound X, the step size is doubled.
318
C The routine is careful never to step beyond SMALL or BIG. If
319
C it hasn't bounded X at SMALL or BIG, QMFINV returns .FALSE.
320
C after setting QLEFT and QHI.
322
C If X is successfully bounded then Algorithm R of the paper
323
C 'Two Efficient Algorithms with Guaranteed Convergence for
324
C Finding a Zero of a Function' by J. C. P. Bus and
325
C T. J. Dekker in ACM Transactions on Mathematical
326
C Software, Volume 1, No. 4 page 330 (DEC. '75) is employed
327
C to find the zero of the function F(X)-Y. This is routine
330
C**********************************************************************
340
STOP '*** EXECUTION FLOWING INTO FLECS PROCEDURES ***'
341
C TO GET-FUNCTION-VALUE
added
removed
Lib/special/cdflib/dinvr.f
