function root_rc ( x, fx, ferr, xerr, q )

c*********************************************************************72
c
cc ROOT_RC solves a single nonlinear equation using reverse communication.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    21 January 2013
c
c  Author:
c
c    Original FORTRAN77 version by Gaston Gonnet.
c    This FORTRAN77 version by John Burkardt.
c
c  Reference:
c     
c    Gaston Gonnet,
c    On the Structure of Zero Finders,
c    BIT Numerical Mathematics,
c    Volume 17, Number 2, June 1977, pages 170-183.
c
c  Parameters:
c
c    Input, double precision X, an estimate for the root.  On the first
c    call, this must be a value chosen by the user.  Thereafter, it may
c    be a value chosen by the user, or the value of ROOT returned on the
c    previous call to the function.
c
c    Input, double precision FX, the value of the function at X.
c
c    Output, double precision FERR, the smallest value of F encountered.
c
c    Output, double precision XERR, the width of the change-in-sign interval,
c    if one was encountered.
c
c    Input/output, double precision Q(9), storage needed by the function.
c    Before the first call, the user must set Q(1) to 0.
c
c    Output, double precision ROOT_RC, an improved estimate for the root.
c
      implicit none

      double precision d
      double precision decr
      double precision ferr
      double precision fx
      integer i
      double precision p
      double precision q(9)
      double precision r
      double precision r8_huge
      double precision r8_sign
      double precision root_rc
      double precision u
      double precision v
      double precision w
      double precision x
      double precision xerr
      double precision z
c
c  If we found an exact zero, there is nothing more to do.
c
      if ( fx .eq. 0.0D+00 ) then
        ferr = 0.0D+00
        xerr = 0.0D+00
        root_rc = x
        return
      end if

      ferr = abs ( fx )
c
c  If this is the first time, initialize, estimate the first root, and exit.
c
      if ( q(1) .eq. 0.0D+00 ) then
        q(1) = fx
        q(2) = x
        do i = 3, 9
          q(i) = 0.0D+00
        end do
        root_rc = x + fx
        xerr = r8_huge ( )
        return
      end if
c
c  This is not the first call.
c
      q(9) = q(9) + 1.0D+00
c
c  Check for too many iterations.
c
      if ( 80.0D+00 .lt. q(9) ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'ROOT_RC - Fatal error!'
        write ( *, '(a,i6)' ) '  Number of iterations = ', int ( q(9) )
        stop
      end if
c
c  Check for a repeated X value.
c
      if ( ( 2.0 .le. q(9) .and. x .eq. q(4) ) .or. 
     &  x .eq. q(2) ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'ROOT_RC - Fatal error!'
        write ( *, '(a,i6)' ) '  Value of X has been input before.'
        stop
      end if
c
c  Push X -> A -> B -> C
c
      do i = 6, 3, -1
        q(i) = q(i-2)
      end do
      q(1) = fx
      q(2) = x
c
c  If we have a change-in-sign interval, store the opposite value.
c
      if ( r8_sign ( q(1) ) .ne. r8_sign ( q(3) ) ) then
        q(7) = q(3)
        q(8) = q(4)
      end if
c
c  Calculate XERR.
c
      if ( q(7) .ne. 0.0D+00 ) then
        xerr = abs ( q(8) - q(2) )
      else
        xerr = r8_huge ( )
      end if
c
c  If more than 30 iterations, and we have change-in-sign interval, bisect.
c
      if ( 30.0D+00 .lt. q(9) .and. q(7) .ne. 0.0D+00 ) then
        root_rc = q(2) + ( q(8) - q(2) ) / 2.0D+00
        return
      end if

      v = ( q(3) - q(1) ) / ( q(4) - q(2) )
c
c  If 3 or more points, try Muller.
c
      if ( q(5) .ne. 0.0D+00 ) then
        u = ( q(5) - q(3) ) / ( q(6) - q(4) )
        w = q(4) - q(2)
        z = ( q(6) - q(2) ) / w
        r = ( z + 1.0D+00 ) * v - u

        if ( r .ne. 0.0D+00 ) then
          p = 2.0D+00 * z * q(1) / r
          d = 2.0D+00 * p / ( w * r ) * ( v - u )
          if ( -1.0D+00 .le. d ) then
            root_rc = q(2) - p / ( 1.0D+00 + sqrt ( 1.0D+00 + d ) )
            if ( q(7) .eq. 0.0D+00 .or.
     &           ( q(2) .lt. root_rc .and. root_rc .lt. q(8) ) .or.
     &           ( q(8) .lt. root_rc .and. root_rc .lt. q(2) ) ) then
              return
            end if
          end if
        end if
      end if
c
c  Try the secant step.
c
      if ( q(1) .ne. q(3) .or. q(7) .eq. 0.0D+00 ) then
        if ( q(1) .eq. q(3) ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'ROOT_RC - Fatal error!'
          write ( *, '(a)' ) '  Cannot apply any method.'
          stop
        end if
        decr = q(1) / v
        if ( abs ( decr ) * 4.6D+18 .lt. abs ( q(2) ) ) then
          decr = 1.74D-18 * abs ( q(2) ) * r8_sign ( decr )
        end if
        root_rc = q(2) - decr
        if ( q(7) .eq. 0.0D+00 .or.
     &      ( q(2) .lt. root_rc .and. root_rc .lt. q(8) ) .or.
     &      ( q(8) .lt. root_rc .and. root_rc .lt. q(2) ) ) then
          return
        end if
      end if
c
c  Apply bisection.
c
      root_rc = q(2) + ( q(8) - q(2) ) / 2.0D+00

      return
      end
      subroutine roots_rc ( n, x, fx, ferr, xnew, q )

c*********************************************************************72
c
cc ROOTS_RC solves a system of nonlinear equations using reverse communication.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    21 January 2013
c
c  Author:
c
c    Original FORTRAN77 version by Gaston Gonnet.
c    This FORTRAN77 version by John Burkardt.
c
c  Reference:
c     
c    Gaston Gonnet,
c    On the Structure of Zero Finders,
c    BIT Numerical Mathematics,
c    Volume 17, Number 2, June 1977, pages 170-183.
c
c  Parameters:
c
c    Input, integer N, the number of equations.
c
c    Input, double precision X(N).  Before the first call, the user should
c    set X to an initial guess or estimate for the root.  Thereafter, the input
c    value of X should be the output value of XNEW from the previous call.
c
c    Input, double precision FX(N), the value of the function at XNEW.
c
c    Output, double precision FERR, the function error, that is, the sum of
c    the absolute values of the most recently computed function vector.
c
c    Output, double precision XNEW(N), a new point at which a function 
c    value is requested.
c
c    Workspace, double precision Q(2*N+2,N+2).  Before the first call 
c    for a given problem, the user must set Q(2*N+1,1) to 0.0.
c
      implicit none

      integer n

      double precision damp
      double precision ferr
      double precision fx(n)
      integer i
      integer info
      integer j
      integer jsma
      integer jsus
      double precision q(2*n+2,n+2)
      double precision r8_huge
      double precision sump
      double precision t
      double precision x(n)
      double precision xnew(n)

      ferr = 0.0D+00
      do i = 1, n
        ferr = ferr + abs ( fx(i) )
      end do
c
c  Initialization if Q(2*N+1,1) = 0.0.
c
      if ( q(2*n+1,1) == 0.0D+00 ) then

        do i = 1, n
          do j = 1, n + 1
            q(i,j) = 0.0D+00
            q(i+1,j) = 0.0D+00
          end do
          q(i,i) = 100.0D+00
          q(i+n,i) = 1.0D+00
        end do
 
        do j = 1, n
          q(2*n+1,j) = r8_huge ( )
        end do

        do j = 1, n
          q(2*n+2,j) = dble ( n )
        end do

        do i = 1, n
          q(i+n,n+1) = x(i)
        end do

        do i = 1, n
          q(i,n+1) = fx(i)
        end do

        q(2*n+1,n+1) = ferr
        q(2*n+2,n+1) = 0.0D+00
        damp = 0.99D+00

      else

        jsus = 1
        do i = 2, n + 1
          if ( dble ( 2 * n ) <= q(2*n+2,i) ) then
            q(2*n+1,i) = r8_huge ( )
          end if
          if ( q(2*n+2,jsus) .lt. ( n + 3 ) / 2 ) then
            jsus = i
          end if
          if ( ( n + 3 ) / 2 .le. q(2*n+2,i) .and. 
     &      q(2*n+1,jsus) .lt. q(2*n+1,i) ) then
            jsus = i
          end if
        end do

        do i = 1, n
          q(i+n,jsus) = x(i)
          q(i,jsus) = fx(i)
        end do

        q(2*n+1,jsus) = ferr
        q(2*n+2,jsus) = 0
        jsma = 1
        damp = 0.0D+00

        do j = 1, n + 1
          if ( r8_huge ( ) / 10.0D+00 .lt. q(2*n+1,j) ) then
            damp = 0.99D+00
          end if
          if ( q(2*n+1,j) .lt. q(2*n+1,jsma) ) then
            jsma = j
          end if
        end do

        if ( jsma .ne. n + 1 ) then
          do i = 1, 2 * n + 2
            t = q(i,jsma)
            q(i,jsma) = q(i,n+1)
            q(i,n+1) = t
          end do
        end if

      end if

      do i = 1, n
        q(i,n+2) = q(i,n+1)
      end do
c
c  Call the linear equation solver, which should not destroy the matrix 
c  in Q(1:N,1:N), and should overwrite the solution into Q(1:N,N+2).
c
      call r8mat_fs ( n, q(1:n,1:n), q(1:n,n+2), info )

      if ( info .ne. 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'ROOTS_RC - Fatal error!'
        write ( *, '(a,i6)' ) 
     &    '  Linear equation solver returns INFO = ', info
        stop
      end if

      sump = 0.0D+00
      do i = 1, n
        sump = sump + q(i,n+2)
      end do

      if ( abs ( 1.0D+00 - sump ) .le. 1.0D-10 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'ROOTS_RC - Fatal error!'
        write ( *, '(a)' ) '  SUMP almost exactly 1.'
        write ( *, '(a,g14.6)' ) '  SUMP = ', sump
        stop
      end if

      do i = 1, n
        xnew(i) = q(i+n,n+1)
        do j = 1, n
          xnew(i) = xnew(i) - q(i+n,j) * q(j,n+2)
        end do
c
c  If system not complete, damp the solution.
c
        xnew(i) = xnew(i) / ( 1.0D+00 - sump ) * ( 1.0D+00 - damp )
     &    + q(i+n,n+1) * damp

      end do

      do j = 1, n + 1
        q(2*n+2,j) = q(2*n+2,j) + 1.0D+00
      end do

      return
      end
      subroutine r8mat_fs ( n, a, b, info )

c*********************************************************************72
c
cc R8MAT_FS factors and solves a system with one right hand side.
c
c  Discussion:
c
c    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
c
c    This routine differs from R8MAT_FSS in two ways:
c    * only one right hand side is allowed;
c    * the input matrix A is not modified.
c
c    This routine uses partial pivoting, but no pivot vector is required.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license. 
c
c  Modified:
c
c    20 January 2013
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, integer N, the order of the matrix.
c    N must be positive.
c
c    Input, double precision A(N,N), the coefficient matrix.
c
c    Input/output, double precision B(N).
c    On input, B is the right hand side of the linear system.
c    On output, B is the solution of the linear system.
c
c    Output, integer INFO, singularity flag.
c    0, no singularity detected.
c    nonzero, the factorization failed on the INFO-th step.
c
      implicit none

      integer n
      integer nb

      double precision a(n,n)
      double precision a2(n,n)
      double precision b(n)
      integer i
      integer info
      integer ipiv
      integer j
      integer jcol
      integer k
      double precision piv
      double precision temp
c
c  Copy the matrix.
c
      do j = 1, n
        do i = 1, n
          a2(i,j) = a(i,j)
        end do
      end do

      info = 0

      do jcol = 1, n
c
c  Find the maximum element in column I.
c
        piv = abs ( a2(jcol,jcol) )
        ipiv = jcol
        do i = jcol + 1, n
          if ( piv .lt. abs ( a2(i,jcol) ) ) then
            piv = abs ( a2(i,jcol) )
            ipiv = i
          end if
        end do

        if ( piv .eq. 0.0D+00 ) then
          info = jcol
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'R8MAT_FS - Fatal error!'
          write ( *, '(a,i8)' ) '  Zero pivot on step ', info
          return
        end if
c
c  Switch rows JCOL and IPIV, and B.
c
        if ( jcol .ne. ipiv ) then

          do j = 1, n
            temp = a2(jcol,j)
            a2(jcol,j) = a2(ipiv,j)
            a2(ipiv,j) = temp
          end do

          temp = b(jcol)
          b(jcol) = b(ipiv)
          b(ipiv) = temp

        end if
c
c  Scale the pivot row.
c
        do j = jcol + 1, n
          a2(jcol,j) = a2(jcol,j) / a2(jcol,jcol)
        end do
        b(jcol) = b(jcol) / a2(jcol,jcol)
        a2(jcol,jcol) = 1.0D+00
c
c  Use the pivot row to eliminate lower entries in that column.
c
        do i = jcol + 1, n
          if ( a2(i,jcol) .ne. 0.0D+00 ) then
            temp = - a2(i,jcol)
            a2(i,jcol) = 0.0D+00
            do j = jcol + 1, n
              a2(i,j) = a2(i,j) + temp * a2(jcol,j)
            end do
            b(i) = b(i) + temp * b(jcol)
          end if
        end do

      end do
c
c  Back solve.
c
      do k = n, 2, -1
        do i = 1, k - 1
          b(i) = b(i) - a2(i,k) * b(k)
        end do
      end do

      return
      end
      function r8_epsilon ( )

c*********************************************************************72
c
cc R8_EPSILON returns the R8 roundoff unit.
c
c  Discussion:
c
c    The roundoff unit is a number R which is a power of 2 with the
c    property that, to the precision of the computer's arithmetic,
c      1 .lt. 1 + R
c    but
c      1 = ( 1 + R / 2 )
c
c    FORTRAN90 provides the superior library routine
c
c      EPSILON ( X )
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    01 September 2012
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Output, double precision R8_EPSILON, the R8 roundoff unit.
c
      implicit none

      double precision r8_epsilon

      r8_epsilon = 2.220446049250313D-016

      return
      end
      function r8_huge ( )

c*********************************************************************72
c
cc R8_HUGE returns a "huge" R8.
c
c  Discussion:
c
c    The value returned by this function is NOT required to be the
c    maximum representable R8.  This value varies from machine to machine,
c    from compiler to compiler, and may cause problems when being printed.
c    We simply want a "very large" but non-infinite number.
c
c    FORTRAN90 provides a built-in routine HUGE ( X ) that
c    can return the maximum representable number of the same datatype
c    as X, if that is what is really desired.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    13 April 2004
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Output, double precision R8_HUGE, a huge number.
c
      implicit none

      double precision r8_huge

      r8_huge = 1.0D+30

      return
      end
      function r8_sign ( x )

c*********************************************************************72
c
cc R8_SIGN returns the sign of an R8.
c
c  Discussion:
c
c    value = -1 if X < 0;
c    value = +1 if X => 0.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    23 August 2008
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, double precision X, the number whose sign is desired.
c
c    Output, double precision R8_SIGN, the sign of X.
c
      implicit none

      double precision r8_sign
      double precision x

      if ( x .lt. 0.0D+00 ) then
        r8_sign = -1.0D+00
      else
        r8_sign = +1.0D+00
      end if

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    12 January 2007
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    None
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 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, month(m), y, h, ':', n, ':', s, '.', mm, ampm

      return
      end