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