1      SUBROUTINE dinvr(status,x,fx,qleft,qhi)
2C**********************************************************************
3C
4C     SUBROUTINE DINVR(STATUS, X, FX, QLEFT, QHI)
5C          Double precision
6C          bounds the zero of the function and invokes zror
7C                    Reverse Communication
8C
9C
10C                              Function
11C
12C
13C     Bounds the    function  and  invokes  ZROR   to perform the   zero
14C     finding.  STINVR  must  have   been  called  before this   routine
15C     in order to set its parameters.
16C
17C
18C                              Arguments
19C
20C
21C     STATUS <--> At the beginning of a zero finding problem, STATUS
22C                 should be set to 0 and INVR invoked.  (The value
23C                 of parameters other than X will be ignored on this cal
24C
25C                 When INVR needs the function evaluated, it will set
26C                 STATUS to 1 and return.  The value of the function
27C                 should be set in FX and INVR again called without
28C                 changing any of its other parameters.
29C
30C                 When INVR has finished without error, it will return
31C                 with STATUS 0.  In that case X is approximately a root
32C                 of F(X).
33C
34C                 If INVR cannot bound the function, it returns status
35C                 -1 and sets QLEFT and QHI.
36C                         INTEGER STATUS
37C
38C     X <-- The value of X at which F(X) is to be evaluated.
39C                         DOUBLE PRECISION X
40C
41C     FX --> The value of F(X) calculated when INVR returns with
42C            STATUS = 1.
43C                         DOUBLE PRECISION FX
44C
45C     QLEFT <-- Defined only if QMFINV returns .FALSE.  In that
46C          case it is .TRUE. If the stepping search terminated
47C          unsuccessfully at SMALL.  If it is .FALSE. the search
48C          terminated unsuccessfully at BIG.
49C                    QLEFT is LOGICAL
50C
51C     QHI <-- Defined only if QMFINV returns .FALSE.  In that
52C          case it is .TRUE. if F(X) .GT. Y at the termination
53C          of the search and .FALSE. if F(X) .LT. Y at the
54C          termination of the search.
55C                    QHI is LOGICAL
56
57C
58C**********************************************************************
59C     Modified by S. Steer INRIA 1998,to replace ASSIGN instruction by
60c     Computed GOTO
61C**********************************************************************
62C     .. Scalar Arguments ..
63      DOUBLE PRECISION fx,x,zabsst,zabsto,zbig,zrelst,zrelto,zsmall,
64     +                 zstpmu
65      INTEGER status
66      LOGICAL qhi,qleft
67C     ..
68C     .. Local Scalars ..
69      DOUBLE PRECISION absstp,abstol,big,fbig,fsmall,relstp,reltol,
70     +                 small,step,stpmul,xhi,xlb,xlo,xsave,xub,yy,zx,zy,
71     +                 zz
72      INTEGER i99999
73      LOGICAL qbdd,qcond,qdum1,qdum2,qincr,qlim,qok,qup
74C     ..
75C     .. External Subroutines ..
76      EXTERNAL dstzr,dzror
77C     ..
78C     .. Intrinsic Functions ..
79      INTRINSIC abs,max,min
80C     ..
81C     .. Statement Functions ..
82      LOGICAL qxmon
83C     ..
84C     .. Save statement ..
85      SAVE
86C     ..
87C     .. Statement Function definitions ..
88      qxmon(zx,zy,zz) = zx .LE. zy .AND. zy .LE. zz
89C     ..
90C     .. Executable Statements ..
91
92      IF (status.GT.0) GO TO 310
93
94      qcond = .NOT. qxmon(small,x,big)
95      IF (qcond) then
96         call basout(io,wte,' SMALL, X, BIG not monotone in INVR')
97         status = -100
98         return
99      endif
100      xsave = x
101C
102C     See that SMALL and BIG bound the zero and set QINCR
103C
104      x = small
105C     GET-FUNCTION-VALUE
106c      ASSIGN 10 TO i99999
107      i99999=1
108      GO TO 300
109
110   10 fsmall = fx
111      x = big
112C     GET-FUNCTION-VALUE
113c      ASSIGN 20 TO i99999
114      i99999=2
115      GO TO 300
116
117   20 fbig = fx
118      qincr = fbig .GT. fsmall
119      IF (.NOT. (qincr)) GO TO 50
120      IF (fsmall.LE.0.0D0) GO TO 30
121      status = -1
122      qleft = .TRUE.
123      qhi = .TRUE.
124      RETURN
125
126   30 IF (fbig.GE.0.0D0) GO TO 40
127      status = -1
128      qleft = .FALSE.
129      qhi = .FALSE.
130      RETURN
131
132   40 GO TO 80
133
134   50 IF (fsmall.GE.0.0D0) GO TO 60
135      status = -1
136      qleft = .TRUE.
137      qhi = .FALSE.
138      RETURN
139
140   60 IF (fbig.LE.0.0D0) GO TO 70
141      status = -1
142      qleft = .FALSE.
143      qhi = .TRUE.
144      RETURN
145
146   70 CONTINUE
147   80 x = xsave
148      step = max(absstp,relstp*abs(x))
149C      YY = F(X) - Y
150C     GET-FUNCTION-VALUE
151c      ASSIGN 90 TO i99999
152      i99999=3
153      GO TO 300
154
155   90 yy = fx
156      IF (.NOT. (yy.EQ.0.0D0)) GO TO 100
157      status = 0
158      qok = .TRUE.
159      RETURN
160
161  100 qup = (qincr .AND. (yy.LT.0.0D0)) .OR.
162     +      (.NOT.qincr .AND. (yy.GT.0.0D0))
163C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
164C
165C     HANDLE CASE IN WHICH WE MUST STEP HIGHER
166C
167C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
168      IF (.NOT. (qup)) GO TO 170
169      xlb = xsave
170      xub = min(xlb+step,big)
171      GO TO 120
172
173  110 IF (qcond) GO TO 150
174C      YY = F(XUB) - Y
175  120 x = xub
176C     GET-FUNCTION-VALUE
177c      ASSIGN 130 TO i99999
178      i99999=4
179      GO TO 300
180
181  130 yy = fx
182      qbdd = (qincr .AND. (yy.GE.0.0D0)) .OR.
183     +       (.NOT.qincr .AND. (yy.LE.0.0D0))
184      qlim = xub .GE. big
185      qcond = qbdd .OR. qlim
186      IF (qcond) GO TO 140
187      step = stpmul*step
188      xlb = xub
189      xub = min(xlb+step,big)
190  140 GO TO 110
191
192  150 IF (.NOT. (qlim.AND..NOT.qbdd)) GO TO 160
193      status = -1
194      qleft = .FALSE.
195      qhi = .NOT. qincr
196      x = big
197      RETURN
198
199  160 GO TO 240
200C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
201C
202C     HANDLE CASE IN WHICH WE MUST STEP LOWER
203C
204C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
205  170 xub = xsave
206      xlb = max(xub-step,small)
207      GO TO 190
208
209  180 IF (qcond) GO TO 220
210C      YY = F(XLB) - Y
211  190 x = xlb
212C     GET-FUNCTION-VALUE
213c      ASSIGN 200 TO i99999
214      i99999=5
215      GO TO 300
216
217  200 yy = fx
218      qbdd = (qincr .AND. (yy.LE.0.0D0)) .OR.
219     +       (.NOT.qincr .AND. (yy.GE.0.0D0))
220      qlim = xlb .LE. small
221      qcond = qbdd .OR. qlim
222      IF (qcond) GO TO 210
223      step = stpmul*step
224      xub = xlb
225      xlb = max(xub-step,small)
226  210 GO TO 180
227
228  220 IF (.NOT. (qlim.AND..NOT.qbdd)) GO TO 230
229      status = -1
230      qleft = .TRUE.
231      qhi = qincr
232      x = small
233      RETURN
234
235  230 CONTINUE
236  240 CALL dstzr(xlb,xub,abstol,reltol)
237C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
238C
239C     IF WE REACH HERE, XLB AND XUB BOUND THE ZERO OF F.
240C
241C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
242      status = 0
243      GO TO 260
244
245  250 IF (.NOT. (status.EQ.1)) GO TO 290
246  260 CALL dzror(status,x,fx,xlo,xhi,qdum1,qdum2)
247      IF (.NOT. (status.EQ.1)) GO TO 280
248C     GET-FUNCTION-VALUE
249c      ASSIGN 270 TO i99999
250      i99999=6
251      GO TO 300
252
253  270 CONTINUE
254  280 GO TO 250
255
256  290 x = xlo
257      status = 0
258      RETURN
259
260      ENTRY dstinv(zsmall,zbig,zabsst,zrelst,zstpmu,zabsto,zrelto)
261C**********************************************************************
262C
263C      SUBROUTINE DSTINV( SMALL, BIG, ABSSTP, RELSTP, STPMUL,
264C     +                   ABSTOL, RELTOL )
265C      Double Precision - SeT INverse finder - Reverse Communication
266C
267C
268C                              Function
269C
270C
271C     Concise Description - Given a monotone function F finds X
272C     such that F(X) = Y.  Uses Reverse communication -- see invr.
273C     This routine sets quantities needed by INVR.
274C
275C          More Precise Description of INVR -
276C
277C     F must be a monotone function, the results of QMFINV are
278C     otherwise undefined.  QINCR must be .TRUE. if F is non-
279C     decreasing and .FALSE. if F is non-increasing.
280C
281C     QMFINV will return .TRUE. if and only if F(SMALL) and
282C     F(BIG) bracket Y, i. e.,
283C          QINCR is .TRUE. and F(SMALL).LE.Y.LE.F(BIG) or
284C          QINCR is .FALSE. and F(BIG).LE.Y.LE.F(SMALL)
285C
286C     if QMFINV returns .TRUE., then the X returned satisfies
287C     the following condition.  let
288C               TOL(X) = MAX(ABSTOL,RELTOL*ABS(X))
289C     then if QINCR is .TRUE.,
290C          F(X-TOL(X)) .LE. Y .LE. F(X+TOL(X))
291C     and if QINCR is .FALSE.
292C          F(X-TOL(X)) .GE. Y .GE. F(X+TOL(X))
293C
294C
295C                              Arguments
296C
297C
298C     SMALL --> The left endpoint of the interval to be
299C          searched for a solution.
300C                    SMALL is DOUBLE PRECISION
301C
302C     BIG --> The right endpoint of the interval to be
303C          searched for a solution.
304C                    BIG is DOUBLE PRECISION
305C
306C     ABSSTP, RELSTP --> The initial step size in the search
307C          is MAX(ABSSTP,RELSTP*ABS(X)). See algorithm.
308C                    ABSSTP is DOUBLE PRECISION
309C                    RELSTP is DOUBLE PRECISION
310C
311C     STPMUL --> When a step doesn't bound the zero, the step
312C                size is multiplied by STPMUL and another step
313C                taken.  A popular value is 2.0
314C                    DOUBLE PRECISION STPMUL
315C
316C     ABSTOL, RELTOL --> Two numbers that determine the accuracy
317C          of the solution.  See function for a precise definition.
318C                    ABSTOL is DOUBLE PRECISION
319C                    RELTOL is DOUBLE PRECISION
320C
321C
322C                              Method
323C
324C
325C     Compares F(X) with Y for the input value of X then uses QINCR
326C     to determine whether to step left or right to bound the
327C     desired x.  the initial step size is
328C          MAX(ABSSTP,RELSTP*ABS(S)) for the input value of X.
329C     Iteratively steps right or left until it bounds X.
330C     At each step which doesn't bound X, the step size is doubled.
331C     The routine is careful never to step beyond SMALL or BIG.  If
332C     it hasn't bounded X at SMALL or BIG, QMFINV returns .FALSE.
333C     after setting QLEFT and QHI.
334C
335C     If X is successfully bounded then Algorithm R of the paper
336C     'Two Efficient Algorithms with Guaranteed Convergence for
337C     Finding a Zero of a Function' by J. C. P. Bus and
338C     T. J. Dekker in ACM Transactions on Mathematical
339C     Software, Volume 1, No. 4 page 330 (DEC. '75) is employed
340C     to find the zero of the function F(X)-Y. This is routine
341C     QRZERO.
342C
343C**********************************************************************
344      small = zsmall
345      big = zbig
346      absstp = zabsst
347      relstp = zrelst
348      stpmul = zstpmu
349      abstol = zabsto
350      reltol = zrelto
351      RETURN
352
353C(jpc)      STOP '*** EXECUTION FLOWING INTO FLECS PROCEDURES ***'
354C     TO GET-FUNCTION-VALUE
355  300 status = 1
356      RETURN
357
358  310 CONTINUE
359      goto(10,20,90,130,200,270) i99999
360c      GO TO i99999
361
362      END
363