hto_root_find2 Module Reference

Functions/Subroutines
real *8 function	hto_zeroin (f, ax, bx, aerr, rerr)

Function/Subroutine Documentation

hto_zeroin()

real*8 function hto_root_find2::hto_zeroin	(		f,
		real*8, intent(in)	ax,
		real*8, intent(in)	bx,
		real*8, intent(in)	aerr,
		real*8, intent(in)	rerr
	)

Definition at line 3428 of file CALLING_cpHTO.f.

*
      USE hto_rootw
* 
! Code converted using TO_F90 by Alan Miller
! Date: 2003-07-14  Time: 12:32:54
 
!-----------------------------------------------------------------------
 
!         FINDING A ZERO OF THE FUNCTION F(X) IN THE INTERVAL (AX,BX)
 
!                       ------------------------
 
!  INPUT...
 
!  F      FUNCTION SUBPROGRAM WHICH EVALUATES F(X) FOR ANY X IN THE
!         CLOSED INTERVAL (AX,BX).  IT IS ASSUMED THAT F IS CONTINUOUS,
!         AND THAT F(AX) AND F(BX) HAVE DIFFERENT SIGNS.
!  AX     LEFT ENDPOINT OF THE INTERVAL
!  BX     RIGHT ENDPOINT OF THE INTERVAL
!  AERR   THE ABSOLUTE ERROR TOLERANCE TO BE SATISFIED
!  RERR   THE RELATIVE ERROR TOLERANCE TO BE SATISFIED
 
!  OUTPUT...
 
!         ABCISSA APPROXIMATING A ZERO OF F IN THE INTERVAL (AX,BX)
 
!-----------------------------------------------------------------------
!  ZEROIN IS A SLIGHTLY MODIFIED TRANSLATION OF THE ALGOL PROCEDURE
!  ZERO GIVEN BY RICHARD BRENT IN ALGORITHMS FOR MINIMIZATION WITHOUT
!  DERIVATIVES, PRENTICE-HALL, INC. (1973).
!-----------------------------------------------------------------------
      
      IMPLICIT NONE
      real*8, INTENT(IN)  :: ax
      real*8, INTENT(IN)  :: bx
      real*8, INTENT(IN)  :: aerr
      real*8, INTENT(IN)  :: rerr
      real*8              :: fn_val
*
* EXTERNAL f
*
      INTERFACE
       FUNCTION f(x) RESULT(fn_val)
       IMPLICIT NONE
       real*8, INTENT(IN)  :: x
       real*8              :: fn_val
       END FUNCTION f
      END INTERFACE
*
      real*8 :: a,b,c,d,e,eps,fa,fb,fc,tol,xm,p,q,r,s,atol,rtol
 
!  COMPUTE EPS, THE RELATIVE MACHINE PRECISION
 
      eps= epsilon(0.0d0)
 
! INITIALIZATION
      
      a= ax
      b= bx
      fa= f(a)
      fb= f(b)
      IF(fa*fb.gt.0.d0) THEN
       inc= 1
      ENDIF
      atol= 0.5d0*aerr
      rtol= max(0.5d0*rerr, 2.0d0*eps)
 
! BEGIN STEP
 
   10 c= a
      fc= fa
      d= b-a
      e= d
   20 IF(abs(fc) < abs(fb)) THEN
         a= b
         b= c
         c= a
         fa= fb
         fb= fc
         fc= fa
      END IF
 
! CONVERGENCE TEST
 
      tol= rtol*max(abs(b),abs(c))+atol
      xm= 0.5d0*(c-b)
      IF(abs(xm) > tol) THEN
       IF(fb /= 0.0d0) THEN
            
! IS BISECTION NECESSARY
    
        IF(abs(e) >= tol) THEN
         IF(abs(fa) > abs(fb)) THEN
        
! IS QUADRATIC INTERPOLATION POSSIBLE
        
          IF(a == c) THEN
          
! LINEAR INTERPOLATION
          
           s= fb/fc
           p= (c-b)*s
           q= 1.0d0-s
          ELSE
          
! INVERSE QUADRATIC INTERPOLATION
          
           q= fa/fc
           r= fb/fc
           s= fb/fa
           p= s*((c-b)*q*(q-r)-(b-a)*(r-1.0d0))
           q= (q-1.0d0)*(r-1.0d0)*(s-1.0d0)
          ENDIF
        
! ADJUST SIGNS
        
          IF(p > 0.0d0) q= -q
          p= abs(p)
        
! IS INTERPOLATION ACCEPTABLE
        
          IF(2.0*p < (3.0*xm*q-abs(tol*q))) THEN
           IF(p < abs(0.5d0*e*q)) THEN
            e= d
            d= p/q
            GO TO 30
           ENDIF
          ENDIF
         ENDIF
        ENDIF
    
! BISECTION
    
        d= xm
        e= d
    
! COMPLETE STEP
    
  30    a= b
        fa= fb
        IF(abs(d) > tol) b= b+d
        IF(abs(d) <= tol) b= b+sign(tol,xm)
        fb= f(b)
        IF(fb*(fc/abs(fc)) > 0.0d0) GO TO 10
        GO TO 20
       ENDIF
      ENDIF
 
! DONE
 
      fn_val= b
      RETURN