zero.f90

Go to the documentation of this file.

subroutine zero(x1, y1, x2, y2, func, err, x, y, izero, ll)

!     ******************************************************************
!     * the zero y=func(x)=0 on the interval (x1,x2) is found.         *
!     * the function func(x) should be provided by an external.        *
!     * upon return izero is the number of iterations which were       *
!     * required to obtain abs(y).le.err.                              *
!     * diagnostic information is printed if l.ne.0.                   *
!     * the argument ll has the double function of communicating the   *
!     * output unit iout=ll/10 and the print switch l=ll-iout*10.      *
!     * iout=0 (l=ll): file is "output".                               *
!     * modified by jan rem to improve convergence 25/08/84.           *
!     ******************************************************************

  use itm_types

  implicit none

  real(r8), intent(inout) :: x1, y1, x2, y2
  real(r8), external :: func
  real(r8), intent(in) :: err
  real(r8), intent(out) :: x, y
  integer(itm_i4), intent(out) :: izero
  integer(itm_i4), intent(in) :: ll

  real(r8) :: x1s, y1s, x2s, y2s, sig, x0, y0, a, b, c, xn, yn
  integer(itm_i4) :: l, iout
  integer(itm_i4) , parameter :: nizero = 50

  iout = ll / 10
  l = ll - iout * 10

  if (iout == 0) iout = 6

  if ((y1 < 0._r8 .and. y2 < 0._r8) .or. (y1 > 0._r8 &
   .and. y2 > 0._r8)) then
     write(iout, 5)
     return
  end if
  if (l /= 0) write(iout, 6)
  x1s = x1
  y1s = y1
  x2s = x2
  y2s = y2
  if (x1 > x2) then
    x1 = x2s
    y1 = y2s
    x2 = x1s
    y2 = y1s
  end if
  sig = 1._r8
  if (y1 >= 0._r8) then
    sig = -1._r8
    y1 = -y1
    y2 = -y2
  end if

!-- begin loop on izero
  izero = 0

  loop_zero : do
    x0 = x1 - (x2 - x1) * y1 / (y2 - y1)
    y0 = sig * func(x0)
    izero = izero + 1
    if (l /= 0) write(iout, 11) izero, x1, y1, x2, y2, x0, y0

    if (abs(y0) <= err) then
      x = x0
      y = y0
      exit
    end if

    if (abs(y0) >= 0.2_r8 * min(-y1, y2)) then
      a = ((y2 - y0) / (x2 - x0) - (y0 - y1) / (x0 - x1)) / (x2 - x1)
      b = (y2 - y1) / (x2 - x1) - a * (x2 + x1)
      c = y0 - a * x0 * x0 - b * x0
      xn = (-b + sqrt(b * b - 4._r8 * a * c)) / (2._r8 * a)
      yn = sig * func(xn)
      izero = izero + 1
      if (l /= 0) write(iout, 12) izero, x1, y1, x2, y2, xn, yn

      if (abs(yn) <= err) then   
        x = xn
        y = yn
        exit
      end if

      if (yn < 0._r8) then
        x1 = xn
        y1 = yn
        if (y0 > 0._r8) then
          x2 = x0
          y2 = y0
        end if
      else
        x2 = xn
        y2 = yn
        if (y0 < 0._r8) then
          x1 = x0
          y1 = y0
        end if
      end if
    else if (y0 > err) then
      x2 = x0
      y2 = y0
    else
      x1 = x0
      y1 = y0
    end if
    if (izero >= nizero) then
      write(iout, 13) nizero 
      x = x0
      y = y0
      exit
    end if
  end do loop_zero
!-- end loop on izero

  x1 = x1s
  y1 = y1s
  x2 = x2s
  y2 = y2s

  return

 5 format(/1x, '***subroutine zero: y1 and y2 violate requirements')
 6 format(/30x, '==== subroutine zero ====')
11 format(1x, 'izero=', i3, 3x, 'x1=', 1pe12.4, 3x, 'y1=', 1pe12.4, &
    3x, 'x2=', 1pe12.4, 3x, 'y2=', 1pe12.4, 3x, 'x0=', 1pe12.4, &
    3x, 'y0=', 1pe12.4)
12 format(1x, 'izero=', i3, 3x, 'x1=', 1pe12.4, 3x, 'y1=', 1pe12.4, &
    3x, 'x2=', 1pe12.4, 3x, 'y2=', 1pe12.4, 3x, 'xn=', 1pe12.4, &
    3x, 'yn=', 1pe12.4)
13 format(/1x, '***subroutine zero: no convergence for x in ', i3, &
    ' steps')

end subroutine zero