program main

!*****************************************************************************80
!
!! zoomin_test() tests zoomin().
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    18 January 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'zoomin_test():'
  write ( *, '(a)' ) '  Fortran90 version'
  write ( *, '(a)' ) '  Test zoomin().'

  call test01 ( )
  call test02 ( )
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'zoomin_test():'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine test01 ( )

!*****************************************************************************80
!
!! test01() runs the tests on a polynomial function.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    18 January 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) a
  real ( kind = rk ) abserr
  real ( kind = rk ) b
  real ( kind = rk ) c
  real ( kind = rk ), external :: func01
  integer ipoly
  integer kmax
  integer mult
  integer nder
  integer nsub

  a = 1.0D+00
  b = 1.5D+00
  c = 4.0D+00
!
!  Set the error tolerance.
!
  abserr = 1.0D-08
!
!  Tell the program how many derivatives are available.
!
  nder = 3
!
!  IPOLY is the order of the polynomial function
!  or -1 if the function is not a polynomial.
!
  ipoly = 3
!
!  KMAX is the maximum number of iterations.
!
  kmax = 30
!
!  MULT is the multiplicity of the root.
!  If not known, set MULT to 1.
!
  mult = 1
!
!  For modified Newton methods, set the number of substeps.
!
  nsub = 3
!
!  Call ZOOMIN.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'test01():'
  write ( *, '(a)' ) '  (Polynomial function F(X))'
  write ( *, '(a)' ) '  Find a root of F(X)=(X+3)*(X+3)*(X-2)=0'
  write ( *, '(a)' ) ' '

  call zoomin_all_test ( a, b, c, abserr, kmax, func01, ipoly, mult, nder, nsub )

  return
end
function func01 ( x, ider )

!*****************************************************************************80
!
!! func01() computes the function value for the first test.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    17 December 1998
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real X, the point at which the evaluation is to take place.
!
!    integer IDER, specifies what is to be evaluated:
!    0, evaluate the function.
!    1, evaluate the first derivative.
!    2, evaluate the second derivative.
!    3, evaluate the third derivative.
!
!  Output:
!
!    real FUNC01, the value of the function or derivative.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) func01
  integer ider
  real ( kind = rk ) x

  if ( ider == 0 ) then
    func01 = ( x + 3.0D+00 )**2 * ( x - 2.0D+00 )
  else if ( ider == 1 ) then
    func01 = ( x + 3.0D+00 ) * ( 3.0D+00 * x - 1.0D+00 )
  else if ( ider == 2 ) then
    func01 = 6.0D+00 * x + 8.0D+00
  else if ( ider == 3 ) then
    func01 = 6.0D+00
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'func01(): Fatal error!'
    write ( *, '(a,i8)' ) '  Derivative of order IDER = ', ider
    write ( *, '(a)' ) '  was requested.'
    stop
  end if

  return
end
subroutine test02 ( )

!*****************************************************************************80
!
!! test02() runs the tests on a nonpolynomial function.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    18 January 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) a
  real ( kind = rk ) abserr
  real ( kind = rk ) b
  real ( kind = rk ) c
  real ( kind = rk ), external :: func02
  integer ipoly
  integer kmax
  integer mult
  integer nder
  integer nsub

  a = 0.9D+00
  b = 0.4D+00
  c = 0.5D+00
!
!  Set the error tolerance.
!
  abserr = 1.0D-08
!
!  Tell the program how many derivatives are available.
!
  nder = 3
!
!  IPOLY is the order of the polynomial function
!  or -1 if the function is not a polynomial.
!
  ipoly = - 1
!
!  KMAX is the maximum number of iterations.
!
  kmax = 60
!
!  MULT is the multiplicity of the root.
!  If not known, just set MULT to 1.
!
  mult = 1
!
!  For modified Newton methods, set the number of substeps.
!
  nsub = 3
!
!  Call ZOOMIN.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'test02():'
  write ( *, '(a)' ) '  (Nonpolynomial function F(X))'
  write ( *, '(a)' ) '  Find a root of F(X) = COS(X) - X'
  write ( *, '(a)' ) ' '

  call zoomin_all_test ( a, b, c, abserr, kmax, func02, ipoly, mult, nder, nsub )

  return
end
function func02 ( x, ider )

!*****************************************************************************80
!
!! func02() computes the function value for the second test.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    17 December 1998
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real X, the point at which the evaluation is to take place.
!
!    integer IDER, specifies what is to be evaluated:
!    0, evaluate the function.
!    1, evaluate the first derivative.
!    2, evaluate the second derivative.
!    3, evaluate the third derivative.
!
!  Output:
!
!    real FUNC02, the value of the function or derivative.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) func02
  integer ider
  real ( kind = rk ) x

  if ( ider == 0 ) then
    func02 = cos ( x ) - x
  else if ( ider == 1 ) then
    func02 = - sin ( x ) - 1.0D+00
  else if ( ider == 2 ) then
    func02 = - cos ( x )
  else if ( ider == 3 ) then
    func02 = sin ( x )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'func02(): Fatal error!'
    write ( *, '(a,i8)' ) '  Derivative of order IDER = ', ider
    write ( *, '(a)' ) '  was requested.'
    stop
  end if

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! timestamp() prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  integer values(8)
  integer y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end
subroutine zoomin_all_test ( a, b, c, abserr, kmax, f, ipoly, mult, nder, nsub )

!*****************************************************************************80
!
!! zoomin_all_test() calls all the zero finders for a given problem.
!
!  Discussion:
!
!    The original code was written by Harold Deiss in BASIC.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    05 August 2005
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real A, B, C, three estimates for a root of the
!    function.  The bisection routines will only be called if the function
!    evaluated at A, B and C includes both positive and negative values.
!
!    real ABSERR, an error tolerance.
!
!    real external F, the name of the routine that
!    evaluates the function or its derivatives, of the form
!      function f ( x, ider )
!
!    integer IPOLY, is 0 if the function is not known to be
!    a polynomial, or is equal to the polynomial degree otherwise.
!
!    integer KMAX, the maximum number of iterations allowed.
!
!    integer MULT, the multiplicity of the root being 
!    sought.  If the multiplicity of the root is not known, set MULT to 1.
!
!    integer NDER, the highest order derivative that can be
!    computed by the user function.  NDER = 0 means only the function
!    value itself is available, for instance.  NDER must be at least 0,
!    and no algorithm needs a derivative higher than 3.
!
!    integer NSUB, the number of substeps to take.  A few
!    algorithms include a "sub-iteration".  For instance, the modified
!    Newton method evaluates the derivative function only before
!    steps 1, NSUB+1, 2*NSUB+1 and so on, and uses that derivative
!    value for NSUB iterations in a row.  This option is useful when
!    the derivative evaluation is expensive.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) a
  real ( kind = rk ) abserr
  real ( kind = rk ) b
  real ( kind = rk ) beta
  real ( kind = rk ) c
  logical change
  character ( len = 3 ) ctemp
  real ( kind = rk ) d
  real ( kind = rk ) e
  real ( kind = rk ) em
  real ( kind = rk ), external :: f
  real ( kind = rk ) fa
  real ( kind = rk ) fb
  real ( kind = rk ) fc
  real ( kind = rk ) ff
  real ( kind = rk ) g
  real ( kind = rk ) h
  integer nder
  integer ierror
  integer ipoly
  integer k
  integer kmax
  integer mult
  character ( len = 80 ) name
  integer nsub
  real ( kind = rk ) rho
  real ( kind = rk ) sa
  real ( kind = rk ) sb
  real ( kind = rk ) sc
  real ( kind = rk ) t
  real ( kind = rk ) x
  real ( kind = rk ) x0
  real ( kind = rk ) x1
  real ( kind = rk ) x2

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'zoomin_all_test():'
  write ( *, '(a)' ) '  A compilation of scalar zero finders,'
  write ( *, '(a)' ) '  based on the work of Joseph Traub.'
  write ( *, '(a)' ) ' '

  fa = f ( a, 0 )
  fb = f ( b, 0 )
  fc = f ( c, 0 )
!
!  Rearrange the data, if necessary, so that the pair (A,FA) has
!  the smallest function value of the three sets of data.
!
  if ( abs ( fb ) < abs ( fa ) ) then
    t = a
    a = b
    b = t
    t = fa
    fa = fb
    fb = t
  end if

  if ( abs ( fc ) < abs ( fa ) ) then
    t = a
    a = c
    c = t
    t = fa
    fa = fc
    fc = t
  end if
!
!  If necessary, switch B and C, so that FB is of opposite sign to FA,
!  if possible.
!
  sa = sign ( 1.0D+00, fa )
  sb = sign ( 1.0D+00, fb )
  sc = sign ( 1.0D+00, fc )

  if ( sa /= sb ) then
    change = .true.
  else if ( sa /= sc ) then
    t = b
    b = c
    c = t
    t = fb
    fb = fc
    fc = t
    change = .true.
  else
    change = .false.
  end if
!
!  Check other input.
!
  if ( mult < 1 ) then
    mult = 1
  end if

  if ( nder < 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'zoomin_all_test(): Fatal error!'
    write ( *, '(a,i4)' ) '  The input quantity NDER = ', nder
    write ( *, '(a)' ) '  But NDER must be at least 0.'
    return
  end if

  if ( nsub < 1 ) then
    nsub = 1
  end if
!
!  Report start up conditions.
!
  write ( *, '(a)' ) '  1 point formulas use:'
  write ( *, '(a,g14.6)' ) '    x1 = ', a
  write ( *, '(a,g14.6)' ) '    fx1=', fa
  write ( * , '(a)' ) '  2 point formulas add:'
  write ( *, '(a,g14.6)' ) '    x2 = ', b
  write ( *, '(a,g14.6)' ) '    fx2=', fb
  write ( * , '(a)' ) '  3 point formulas add:'
  write ( *, '(a,g14.6)' ) '    x3 = ', c
  write ( *, '(a,g14.6)' ) '    fx3=', fc
  write ( * , '(a)' ) ' '
  write ( * , '(a,i4)' ) '  User estimated multiplicity = ', mult
  if ( 0 <= ipoly ) then
    write ( * , '(a,i4)' ) '  Polynomial degree = ', ipoly
  else
    write ( *, '(a)' ) '  The function is not known to be polynomial.'
  end if
  write ( * , '(a,i4)' ) '  Highest derivative supplied = ', nder
  write ( * , '(a,g14.6)' ) '  Error tolerance = ', abserr
  write ( * , '(a,i4)' ) '  Maximum number of steps = ', kmax
  write ( * , '(a,i4)' ) '  Newton method substep parameter  = ', nsub
  write ( * , '(a)' ) ' '
  write ( * , '(a)' ) ' '
  write ( * , '(a)' ) &
    '  Technique                           Root        Steps' &
     // '   Error     Multiplicity'

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  1. One point iteration functions.'
  write ( *, '(a)' ) ' '
!
!  Newton method.
!
  if ( 1 <= nder ) then

    name = 'Newton'
    x = a
    call newton ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Steffenson's method.
!
  name = 'Steffenson'
  x = a
  call steffenson ( x, abserr, kmax, f, ierror, k )

  if ( ierror /= 0 ) then
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
  else
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
  end if
!
!  Stirling's method.
!
  if ( 1 <= nder ) then

    name = 'Stirling'
    x = a
    call stirling ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Midpoint method.
!
  if ( 1 <= nder ) then

    name = 'midpoint'
    x = a
    call midpoint ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Traub-Ostrowski method.
!
  name = 'Traub-Ostrowski'
  x = a
  call traub_ostrowski ( x, abserr, kmax, f, ierror, k )

  if ( ierror /= 0 ) then
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
  else
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
  end if
!
!  Chebyshev method.
!
  if ( 2 <= nder ) then

    name = 'Chebyshev'
    x = a
    call chebyshev ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Halley Super method.
!
  if ( 2 <= nder ) then

    name = 'Halley Super'
    x = a
    call halley_super ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Whittaker convex acceleration method.
!
  if ( 2 <= nder ) then

    name = 'Whittaker'
    x = a
    call whittaker ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Whittaker double convex acceleration method.
!
  if ( 2 <= nder ) then

    name = 'Whittaker2'
    x = a
    call whittaker2 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  E3.
!
  if ( 2 <= nder ) then

    name = 'E3'
    x = a
    call e3 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  E4.
!
  if ( 3 <= nder ) then

    name = 'E4'
    x = a
    call e4 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Halley's method.
!
  if ( 2 <= nder ) then

    name = 'Halley'
    x = a
    call halley1 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  PSI(2,1).
!
  if ( 3 <= nder ) then

    name = 'Psi(2,1)'
    x = a
    call psi_21 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  PSI(1,2).
!
  if ( 3 <= nder ) then

    name = 'Psi(1,2)'
    x = a
    call psi_i2 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Capital PHI(0,3).
!
  if ( 2 <= nder ) then

    name = 'Capital Phi(0,3)'
    x = a
    call cap_phi_03 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Reduced capital PHI(0,4).
!
  if ( 3 <= nder ) then

    name = 'Reduced Capital Phi(0,4)'
    x = a
    call red_cap_phi_04 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Functions from Hansen and Patrick.
!
!  The Ostrowski square root formula is the Hansen family member with BETA = 0.
!
    if ( 2 <= nder ) then
      name = 'Ostrowski square root'
      x = a
      call ostrowski_sqrt ( x, abserr, kmax, f, ierror, k )

      if ( ierror /= 0 ) then
        write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
      else
        write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
      end if

    end if
!     end if
!
!  The Euler method.
!
  if ( 2 <= nder ) then

    name = 'Euler'
    x = a
    call euler ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  The Laguerre method.
!
  if ( 2 <= nder .and. 2 <= ipoly ) then

    name = 'Laguerre'
    x0 = a
    call laguerre ( x0, ipoly, abserr, kmax, f, x, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if

  write ( * , '(a)' ) ' '
  write ( * , '(a)' ) '  2.  One point iteration functions with memory:'
  write ( * , '(a)' ) ' '
!
!  Secant algorithm.
!
  name = 'Secant'
  x = a
  x1 = b
  call secant ( x, x1, abserr, kmax, f, ierror, k )

  if ( ierror /= 0 ) then
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
  else
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
  end if
!
!  Extended secant algorithm.
!
  name = 'Extended secant'
  x = a
  x1 = b
  x2 = c
  call secant_extended ( x, x1, x2, abserr, kmax, f, ierror, k )

  if ( ierror /= 0 ) then
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
  else
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
  end if
!
!  Capital Phi(2,1).
!
  name = 'Capital Phi(2,1)'
  x = a
  x1 = b
  x2 = c
  call cap_phi_21 ( x, x1, x2, abserr, kmax, f, ierror, k )

  if ( ierror /= 0 ) then
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
  else
    write ( *, '(2x,a30,g15.6, i8 )' ) name, x, k
  end if
!
!  Muller's algorithm.
!
  name = 'R8_Muller'
  x = a
  x1 = b
  x2 = c
  call r8_muller ( x, x1, x2, abserr, kmax, f, ierror, k )

  if ( ierror /= 0 ) then
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
  else
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
  end if
!
!  Perp E(2,1).
!
  if ( 1 <= nder ) then

    name = 'Perp E(2,1)'
    x = a
    x1 = b
    x2 = c
    call perp_e_21 ( x, x1, x2, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Star E(2,1).
!
  name = 'Star E(2,1)'
  x = a
  x1 = b
  x2 = c
  call star_e21 ( x, x1, x2, abserr, kmax, f, ierror, k )

  if ( ierror /= 0 ) then
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
  else
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
  end if
!
!  Finite difference Halley's method.
!
  name = 'Finite difference Halley'
  x = a
  x1 = b
  x2 = c
  call halley2 ( x, x1, x2, abserr, kmax, f, ierror, k )

  if ( ierror /= 0 ) then
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
  else
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
  end if
!
!  Phi(1,2).
!
  if ( 1 <= nder ) then

    name = 'Phi(1,2)'
    x = a
    x1 = b
    call phi_12 ( x, x1, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Perp E(1,2).
!
  if ( 1 <= nder ) then

    name = 'Perp E(1,2)'
    x = a
    x1 = b
    call perp_e_12 ( x, x1, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Star E(1,2).
!
  if ( 1 <= nder ) then

    name = 'Star E(1,2)'
    x = a
    x1 = b
    call star_e12 ( x, x1, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Dagger E(1,2).
!
  if ( 1 <= nder ) then

    name = 'Dagger E(1,2)'
    x = a
    x1 = b
    call dagger_e12 ( x, x1, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '  3.  Multiple root methods, multiplicity given.'
  write ( *, '(a)' ) ' '
!
!  Script E2.
!
  if ( 1 <= nder ) then

    name = 'Script E2'
    x = a
    call script_e2 ( x, mult, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Script E3.
!
  if ( 2 <= nder ) then

    name = 'Script E3'
    x = a
    call script_e3 ( x, mult, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Script E4.
!
  if ( 3 <= nder ) then

    name = 'Script E4'
    x = a
    call script_e4 ( x, mult, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Star E - 1,1(f) function in Traub, page 235.
!
  name = 'Star E 1,1(f)'
  x = a
  x1 = b
  call te11f ( x, x1, mult, abserr, kmax, f, ierror, k )

  if ( ierror /= 0 ) then
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
  else
    write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
  end if
!
!  Routines which determine the root and multiplicity
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  4.  Multiple root methods,'
  write ( *, '(a)' ) '      multiplicity not given.'
  write ( *, '(a)' ) ' '
!
!  E2(U).
!
  if ( 2 <= nder ) then

    name = 'E2(U)'
    x = a
    call te2u ( x, em, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, 5x, g15.6 )' ) name, x, k, em
    end if

  end if
!
!  Phi - 1,1(u) function in Traub, page 235.
!
  if ( 1 <= nder ) then

    name = 'Phi 1,1(U)'
    x = a
    call tphi1u ( x, em, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, 5x, g15.6 )' ) name, x, k, em
    end if

  end if
!
!  Third function in Traub, page 235 bottom.
!
  if ( 1 <= nder ) then

    name = 'Traub third'
    x = a
    call tthip ( x, em, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, 5x, g15.6 )' ) name, x, k, em
    end if

  end if
!
!  Van de Vel iteration.
!
  if ( 1 <= nder ) then
    name = 'Van de Vel'
    x = a
    call vandevel ( x, em, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, 5x, g15.6 )' ) name, x, k, em
    end if

  end if
!
!  Improved Van de Vel iteration.
!
  if ( 1 <= nder ) then
    name = 'Van de Vel improved'
    x = a
    call vandevel_improved ( x, em, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, 5x, g15.6 )' ) name, x, k, em
    end if

  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  5.  Multipoint iteration functions.'
  write ( *, '(a)' ) ' '
!
!  First function in Traub, page 236.
!
  if ( 1 <= nder ) then

    name = 'Traub first'
    x = a
    call traub1 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  A family of iteration functions including
!  forms 2, 12, and 13 of Traub, pages 236 - 237.
!
  if ( 1 <= nder ) then

    name = 'Traub second'
    e = 0.5D+00
    d = 1.0D+00
    x = a
    call t_family1 ( x, e, d, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if

  if ( 1 <= nder ) then

    name = 'Traub twelfth'
    e = 0.25D+00
    d = 2.0D+00 / 3.0D+00
    x = a
    call t_family1 ( x, e, d, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if

  if ( 1 <= nder ) then

    name = 'Traub thirteenth'
    e = 5.0D+00 / 12.0D+00
    d = 6.0D+00 / 7.0D+00
    x = a
    call t_family1 ( x, e, d, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Modified Newton method.
!
  if ( 1 <= nder .and. 1 <= nsub ) then

    write ( name, '(a,i3)' ) 'Modified Newton, NSUB=', nsub
    x = a
    call newton_mod ( x, nsub, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Fourth function in Traub.
!
  if ( 1 <= nder .and. 1 <= nsub ) then

    name = 'Traub fourth, NSUB='
    write (ctemp,'(i3)') nsub
    name(20:22) = ctemp(1:3)
    x = a
    call traub4 ( x, nsub, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Newton - secant.
!
  if ( 1 <= nder ) then

    name = 'Newton - secant'
    x = a
    call newton_secant ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  A family of iteration functions including
!  forms 6 and 7 of Traub, pages 236 - 237.
!
  if ( 1 <= nder ) then

    name = 'Traub sixth'
    e = 2.0D+00
    ff = 3.0D+00
    g = 1.0D+00
    h = 1.0D+00
    x = a
    call t_family2 ( x, e, ff, g, h, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if

  if ( 1 <= nder ) then

    name = 'Traub seventh'
    e = 4.0D+00
    ff = 7.0D+00
    g = 3.0D+00
    h = 2.0D+00 / 3.0D+00
    x = a
    call t_family2 ( x, e, ff, g, h, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Eighth function in Traub, page 237.
!
  if ( 1 <= nder ) then

    name = 'Traub eighth'
    x = a
    call traub8 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Ninth function in Traub, page 237.
!
  if ( 1 <= nder ) then

    name = 'Traub ninth'
    x = a
    call traub9 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Type 1 family of iteration functions including
!  forms 10 and 11 of Traub, pages 237.
!
  if ( 1 <= nder ) then

    name = 'Traub type 1, form 10'
    rho = ( 1.0D+00 - sqrt ( 5.0D+00 ) ) / 2.0D+00
    e = 0.0D+00
    x = a
    call tt1f ( x, e, rho, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if

  if ( 1 <= nder ) then

    name = 'Traub type 1, form 11'
    e = 1.0D+00
    x = a
    call tt1f ( x, e, rho, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Fourteenth function in Traub, page 237
!
  if ( 1 <= nder ) then

    name = 'Traub fourteenth'
    x = a
    call t14 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Fifteenth function in Traub, page 238.
!
  if ( 1 <= nder ) then

    name = 'Traub fifteenth'
    x = a
    call t15 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Sixteenth function in Traub, page 238.
!
  if ( 1 <= nder ) then

    name = 'Traub sixteenth'
    x = a
    call t16 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  A family of fourth order methods by Richard King.
!
  if ( 1 <= nder ) then

    name = 'King, BETA=0'
    beta = 0.0D+00
    x = a
    call king ( x, beta, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if

  if ( 1 <= nder ) then

    name = 'King, BETA=1'
    beta = 1.0D+00
    x = a
    call king ( x, beta, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if

  if ( 1 <= nder ) then

    name = 'King, BETA=2'
    beta = 2.0D+00
    x = a
    call king ( x, beta, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Jarratt iterative technique.
!
  if ( 1 <= nder ) then

    name = 'Jarratt'
    x = a
    call jarratt ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Jarratt inverse-free iterative technique.
!
  if ( 1 <= nder ) then

    name = 'Jarratt inverse-free'
    x = a
    call jarratt2 ( x, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if
!
!  Bisection methods
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  6. Bisection methods'
  write ( *, '(a)' ) ' '
!
!  Bisection methods may only be used if there is a change of sign.
!
  if ( .not. change ) then

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'zoomin_all_test():'
    write ( *, '(a)' ) '  Bisection methods will not be called since'
    write ( *, '(a)' ) '  there is no change of sign interval.'

  else
!
!  Bisection.
!
    name = 'Bisection'
    x = a
    x1 = b
    call bisect ( x, x1, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if
!
!  Regula falsi.
!
    name = 'Regula falsi'
    x = a
    x1 = b
    call regula ( x, x1, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if
!
!  Brent.
!
    name = 'Brent'
    x = a
    x1 = b
    call brent ( x, x1, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if
!
!  Bisection + secant.
!
    name = 'Bisection + secant'
    x = a
    x1 = b
    call rhein1 ( x, x1, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if
!
!  Bisection + secant + inverse quadratic.
!
    name = 'Bisection + secant + inv quad'
    x = a
    x1 = b
    call rhein2 ( x, x1, abserr, kmax, f, ierror, k )

    if ( ierror /= 0 ) then
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k, ierror
    else
      write ( *, '(2x,a30,g15.6, i8, i5)' ) name, x, k
    end if

  end if

  return
end