c To get r1mach, mail netlib c send r1mach from core real function seroin(ax,bx,f,tol) real ax,bx,f,tol c c a zero of the function f(x) is computed in the interval ax,bx . c c input.. c c ax left endpoint of initial interval c bx right endpoint of initial interval c f function subprogram which evaluates f(x) for any x in c the interval ax,bx c tol desired length of the interval of uncertainty of the c final result (.ge.0.) c c output.. c c zeroin abscissa approximating a zero of f in the interval ax,bx c c it is assumed that f(ax) and f(bx) have opposite signs c this is checked, and an error message is printed if this is not c satisfied. seroin returns a zero x in the given interval c ax,bx to within a tolerance 4*macheps*abs(x)+tol, where macheps is c the relative machine precision defined as the smallest representable c number such that 1.+macheps .gt. 1. c this function subprogram is a slightly modified translation of c the algol 60 procedure zero given in richard brent, algorithms for c minimization without derivatives, prentice-hall, inc. (1973). c real a,b,c,d,e,eps,fa,fb,fc,tol1,xm,p,q,r,s real abs, r1mach 10 eps = r1mach(4) tol1 = eps+1.0e0 c a=ax b=bx fa=f(a) fb=f(b) c check that f(ax) and f(bx) have different signs if (fa .eq.0.0e0 .or. fb .eq. 0.0e0) go to 20 if (fa * (fb/ abs(fb)) .le. 0.0e0) go to 20 write(6,2500) 2500 format(1x,'f(ax) and f(bx) do not have different signs,', 1 ' seroin is aborting') return 20 c=a fc=fa d=b-a e=d 30 if ( abs(fc).ge. abs(fb)) go to 40 a=b b=c c=a fa=fb fb=fc fc=fa 40 tol1=2.0e0*eps* abs(b)+0.5e0*tol xm = 0.5e0*(c-b) if (( abs(xm).le.tol1).or.(fb.eq.0.0e0)) go to 150 c c see if a bisection is forced c if (( abs(e).ge.tol1).and.( abs(fa).gt. abs(fb))) go to 50 d=xm e=d go to 110 50 s=fb/fa if (a.ne.c) go to 60 c c linear interpolation c p=2.0e0*xm*s q=1.0e0-s go to 70 c c inverse quadratic interpolation c 60 q=fa/fc r=fb/fc p=s*(2.0e0*xm*q*(q-r)-(b-a)*(r-1.0e0)) q=(q-1.0e0)*(r-1.0e0)*(s-1.0e0) 70 if (p.le.0.0e0) go to 80 q=-q go to 90 80 p=-p 90 s=e e=d if (((2.0e0*p).ge.(3.0e0*xm*q- abs(tol1*q))).or.(p.ge. * abs(0.5e0*s*q))) go to 100 d=p/q go to 110 100 d=xm e=d 110 a=b fa=fb if ( abs(d).le.tol1) go to 120 b=b+d go to 140 120 if (xm.le.0.0e0) go to 130 b=b+tol1 go to 140 130 b=b-tol1 140 fb=f(b) if ((fb*(fc/ abs(fc))).gt.0.0e0) go to 20 go to 30 150 seroin=b return end