zeroin.f

Go to the documentation of this file.

c
c
c
c     =================================================
      function zeroin(ax,bx,f,tol)        
c     =================================================
      implicit double precision (a-h,o-z)
      external f
c  
c      a zero of the function  f(x)  is computed in the interval ax,bx .
c      (Standard routine from netlib)
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.d0)  
c  
c  
c  output..   
c  
c  zeroin abcissa approximating a zero of  f  in the interval ax,bx  
c  
c  
c      it is assumed  that   f(ax)   and   f(bx)   have  opposite  signs
c  without  a  check.  zeroin  returns a zero  x  in the given interval 
c  ax,bx  to within a tolerance  4*macheps*dabs(x) + tol, where macheps 
c  is the relative machine precision.          
c
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  
c  
c  
c  compute eps, the relative machine precision 
c  
      eps = 1.d0          
   10 eps = eps/2.d0      
      tol1 = 1.d0 + eps   
      if (tol1 .gt. 1.d0) go to 10   
c  
c initialization         
c  
      a = ax  
      b = bx  
      fa = f(a)          
      fb = f(b)          
c  
c begin step  
c  
   20 c = a   
      fc = fa 
      d = b - a          
      e = d   
   30 if (dabs(fc) .ge. dabs(fb)) go to 40       
      a = b   
      b = c   
      c = a   
      fa = fb 
      fb = fc 
      fc = fa 
c  
c convergence test       
c  
   40 tol1 = 2.d0*eps*dabs(b) + 0.5*tol          
      xm = .5*(c - b)    
      if (dabs(xm) .le. tol1) go to 90          
      if (fb .eq. 0.d0) go to 90     
c  
c is bisection necessary 
c  
      if (dabs(e) .lt. tol1) go to 70
      if (dabs(fa) .le. dabs(fb)) go to 70       
c  
c is quadratic interpolation possible          
c  
      if (a .ne. c) go to 50        
c  
c linear interpolation   
c  
      s = fb/fa          
      p = 2.d0*xm*s       
      q = 1.d0 - s        
      go to 60
c  
c inverse quadratic interpolation   
c  
   50 q = fa/fc          
      r = fb/fc          
      s = fb/fa          
      p = s*(2.d0*xm*q*(q - r) - (b - a)*(r - 1.d0))        
      q = (q - 1.d0)*(r - 1.d0)*(s - 1.d0)        
c  
c adjust signs
c  
   60 if (p .gt. 0.d0) q = -q        
      p = dabs(p)         
c  
c is interpolation acceptable       
c  
      if ((2.d0*p) .ge. (3.d0*xm*q - dabs(tol1*q))) go to 70 
      if (p .ge. dabs(0.5*e*q)) go to 70        
      e = d   
      d = p/q 
      go to 80
c  
c bisection   
c  
   70 d = xm  
      e = d   
c  
c complete step          
c  
   80 a = b   
      fa = fb 
      if (dabs(d) .gt. tol1) b = b + d          
      if (dabs(d) .le. tol1) b = b + dsign(tol1, xm)        
      fb = f(b)          
      if ((fb*(fc/dabs(fc))) .gt. 0.d0) go to 20 
      go to 30
c  
c done        
c  
   90 zeroin = b         
      return  
      end