CHALLENGE using root_fortran to find root of a complex function

32

Thanks for all the comments. Although I am new to Fortran, I have been stuck for the last 10 days trying to find root of a function.

I have been trying to try to link root fortran to my work to find the root a function. I use other libraries(such as lapack and also solved root problems with root_fortran on other codes). but I dont know why the root fortran does not work with my code, I also tried debugging. Therefore, I have resolved to writing a similar function to the uniroot function in R. Here is the link to the zeroin https://www.netlib.org/fmm/zeroin.f. but again I guess I am missing a step. could someone please help me look at this in relation to the zeroin function and provide a feedback? I will be very grateful for help

module root_finder
use iso_fortran_env, only: wp => real64
implicit none
contains

function zeroin(ax,bx,f,tol) result(root)
! Inputs
real(wp), intent(in) :: ax, bx
interface
real(wp) function f(x)
import :: wp
real(wp), intent(in) :: x
end function f
end interface
real( wp), intent(in) :: tol
real( wp) :: root
integer :: i, maxiter = 1000

    ! Local variables
    real(wp) :: a,b,c,d,e,eps,fa,fb,fc
    real(wp) :: p,q,r,s,tol1,xm
    
    ! Initialize
    eps = epsilon(1.0_wp)
    eps = eps/2.0_wp
    tol1 = 1.0_wp + eps

    ! Set brackets
    a = ax
    b = bx
    fa = f(a)
    fb = f(b)

    ! Bisection loop
    do i=1,maxiter
      c = a 
      fc = fa
      d = b - a
      e = d

      ! Update brackets
      if (abs(fc) >= abs(fb)) then
         a = b
         b = c
         c = a
         fa = fb
         fb = fc
         fc = fa
      end if

      ! Convergence test
      tol1 = 2.0_wp*eps*abs(b) + 0.5_wp*tol
      xm = 0.5_wp*(c - b)
      if (abs(xm) <= tol1 .or. fb == 0.0_wp) exit

      ! Try interpolation
      if (abs(e) < tol1 .and. abs(fa) < abs(fb)) then
        d = xm
        e = d
        else if (a /= c) then
           q = fa/fc
           r = fb/fc
           s = fb/fa
           p = s*(2.0_wp*xm*q*(q - r) - (b - a)*(r - 1.0_wp)) 
           q = (q - 1.0_wp)*(r - 1.0_wp)*(s - 1.0_wp)
           
        else 
           s = fb/fa
           p = 2.0_wp*xm*s
           q = 1.0_wp - s
         ! Use bisection
         
      end if

      ! Adjust signs and complete step
      if (p > 0.0_wp) q = -q
      p = abs(p)

      if ( (2.0_wp*p) >= (3.0_wp*xm*q - abs(tol1*q)) .and. &
           p >= abs(0.5_wp*e*q) ) then
         ! Interpolation acceptable
         e = d  
         d = p/q     
      else
         ! Use bisection
         d = xm
         e = d
      end if

      ! Update brackets
      a = b
      fa = fb
      if (abs(d) > tol1) then
         b = b + d
      else
         b = b + sign(tol1,xm)
      end if
      fb = f(b)

      if (fb*(fc/abs(fc)) > 0.0_wp) cycle

    end do
    
    ! Return result
    root = b
    
  end function
end module