global_bis_sec Subroutine

    SUBROUTINE global_bis_sec   ( dx      , error   , fnctn   , maxerr  ,      &
                                  maxsteps, numroots, roots   , steps   ,      &
                                  verbose , xmax    , xmin     )
      !! A subroutine which finds the values of all roots of a
      !! function (or as many as will fit into the provided arrays)
      !! within a given range of values. This subroutine uses the
      !! hybrid bisection-secant root-finding algorithm.
        IMPLICIT NONE

        interface
            pure function fnctn ( x        )
              !! The Fortran function for which the root will be
              !! found. Must take only one real argument, return a
              !! real argument.
                real(8), intent(in) ::  x
                real(8)             ::  fnctn
            end function fnctn
        end interface
        
        ! Input and output variables:
        integer, intent(in)                 ::  maxsteps
          !! An integer value for the maximum number of iterations
          !! applied to the algorithm before terminating.
        logical, intent(in)                 ::  verbose
          !! A logical variable which specifies whether to print
          !! progress to stdout as brackets found and at each
          !! iteration as root found. Also says whether to print a
          !! warning if 'dx' set to 'dxmin' during 'globrack' routine
          !! and if maximum number of iterations reached while finding
          !! root.
        real(8), intent(in)                 ::  maxerr
          !! A real value which specifies the the maximum allowable
          !! error in the computed root. The subroutine terminates
          !! once the error passes below this level.
        real(8), intent(in)                 ::  xmax
          !! The upper limit of the range on which the subroutine will
          !! search for roots and brackets.
        real(8), intent(in)                 ::  xmin
          !! The lower limit of the range on which the subroutine will
          !! search for roots and brackets.
        real(8), intent(inout)              ::  dx
          !! The initial size of increment to use when examining
          !! function.  Minimum interval will be 0.01 of this.
        integer, intent(out)                ::  numroots
          !! An integer value which will return the number of roots
          !! for which brackets were found. A negative number
          !! indicates that an error occurred.
        integer, intent(out), dimension(:)  ::  steps
          !! An integer array in which the number of iterations needed
          !! before error falls below maxerr for each root is stored
          !! and returned.
        real(8), intent(out), dimension(:)  ::  error
          !! A real array in which an estimate of the error in each
          !! computed root will be stored and returned.
        real(8), intent(out), dimension(:)  ::  roots
          !! A real array in which the computed values of each root
          !! will be returned.
        
        ! Other variables:
        INTEGER                                 ::  counter = 0,       &
                                                    maxroots = 0
        REAL(8), DIMENSION(:,:), ALLOCATABLE    ::  brackets
    !------------------------------------------------------------------!
    
        ! Maximum number of roots is given by the size of the smallest
        ! of the arrays 'error', 'roots', and 'steps' (of course, if
        ! it is written properly then the calling program should have
        ! made these all the same size)
        maxroots = MIN(SIZE(error),SIZE(roots),SIZE(steps))
        ALLOCATE(brackets(2,maxroots))
        brackets = 0.d0
        
        ! Get the sets of bracket pairs for the roots of 'fnctn', using
        ! the 'globrack' subroutine 
        CALL global_brackets(brackets, dx, fnctn, numroots, verbose,   &
                             xmax, xmin)
        
        ! For each bracket pair, find the enclosed root using the 
        ! 'biscnt' subroutine
        DO counter = 1,numroots
            CALL bis_secant(error(counter), fnctn, maxerr, maxsteps,   &
                            steps(counter), roots(counter),            &
                            brackets(1,counter), brackets(2,counter))
        END DO

        RETURN
    END SUBROUTINE global_bis_sec