subroutine nelmin ( fn, n, start, xmin, ynewlo, reqmin, step, konvge, kcount, &
  icount, numres, ifault )

!*****************************************************************************80
!
!! NELMIN minimizes a function using the Nelder-Mead algorithm.
!
!  Discussion:
!
!    This routine seeks the minimum value of a user-specified function.
!
!    Simplex function minimisation procedure due to Nelder+Mead(1965),
!    as implemented by O'Neill(1971, Appl.Statist. 20, 338-45), with
!    subsequent comments by Chambers+Ertel(1974, 23, 250-1), Benyon(1976,
!    25, 97) and Hill(1978, 27, 380-2)
!
!    The function to be minimized must be defined by a function of
!    the form
!
!      function fn ( x, f )
!      real ( kind = 8 ) fn
!      real ( kind = 8 ) x(*)
!
!    and the name of this subroutine must be declared EXTERNAL in the
!    calling routine and passed as the argument FN.
!
!    This routine does not include a termination test using the
!    fitting of a quadratic surface.
!
!  Modified:
!
!    27 February 2008
!
!  Author:
!
!    FORTRAN77 version by R ONeill
!    FORTRAN90 version by John Burkardt
!
!  Reference:
!
!    John Nelder, Roger Mead,
!    A simplex method for function minimization,
!    Computer Journal,
!    Volume 7, 1965, pages 308-313.
!
!    R ONeill,
!    Algorithm AS 47:
!    Function Minimization Using a Simplex Procedure,
!    Applied Statistics,
!    Volume 20, Number 3, 1971, pages 338-345.
!
!  Parameters:
!
!    Input, external FN, the name of the function which evaluates
!    the function to be minimized.
!
!    Input, integer ( kind = 4 ) N, the number of variables.
!
!    Input/output, real ( kind = 8 ) START(N).  On input, a starting point
!    for the iteration.  On output, this data may have been overwritten.
!
!    Output, real ( kind = 8 ) XMIN(N), the coordinates of the point which
!    is estimated to minimize the function.
!
!    Output, real ( kind = 8 ) YNEWLO, the minimum value of the function.
!
!    Input, real ( kind = 8 ) REQMIN, the terminating limit for the variance
!    of function values.
!
!    Input, real ( kind = 8 ) STEP(N), determines the size and shape of the
!    initial simplex.  The relative magnitudes of its elements should reflect
!    the units of the variables.
!
!    Input, integer ( kind = 4 ) KONVGE, the convergence check is carried out 
!    every KONVGE iterations.
!
!    Input, integer ( kind = 4 ) KCOUNT, the maximum number of function 
!    evaluations.
!
!    Output, integer ( kind = 4 ) ICOUNT, the number of function evaluations 
!    used.
!
!    Output, integer ( kind = 4 ) NUMRES, the number of restarts.
!
!    Output, integer ( kind = 4 ) IFAULT, error indicator.
!    0, no errors detected.
!    1, REQMIN, N, or KONVGE has an illegal value.
!    2, iteration terminated because KCOUNT was exceeded without convergence.
!
  implicit none

  integer ( kind = 4 ) n

  real    ( kind = 8 ), parameter :: ccoeff = 0.5D+00
  real    ( kind = 8 ) del
  real    ( kind = 8 ) dn
  real    ( kind = 8 ) dnn
  real    ( kind = 8 ), parameter :: ecoeff = 2.0D+00
  real    ( kind = 8 ), parameter :: eps = 0.001D+00
  real    ( kind = 8 ), external :: fn
  integer ( kind = 4 ) i
  integer ( kind = 4 ) icount
  integer ( kind = 4 ) ifault
  integer ( kind = 4 ) ihi
  integer ( kind = 4 ) ilo
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jcount
  integer ( kind = 4 ) kcount
  integer ( kind = 4 ) konvge
  integer ( kind = 4 ) l
  integer ( kind = 4 ) nn
  integer ( kind = 4 ) numres
  real    ( kind = 8 ) p(n,n+1)
  real    ( kind = 8 ) p2star(n)
  real    ( kind = 8 ) pbar(n)
  real    ( kind = 8 ) pstar(n)
  real    ( kind = 8 ), parameter :: rcoeff = 1.0D+00
  real    ( kind = 8 ) reqmin
  real    ( kind = 8 ) rq
  real    ( kind = 8 ) start(n)
  real    ( kind = 8 ) step(n)
  real    ( kind = 8 ) x
  real    ( kind = 8 ) xmin(n)
  real    ( kind = 8 ) y(n+1)
  real    ( kind = 8 ) y2star
  real    ( kind = 8 ) ylo
  real    ( kind = 8 ) ynewlo
  real    ( kind = 8 ) ystar
  real    ( kind = 8 ) z
!
!  Check the input parameters.
!
  if ( reqmin <= 0.0D+00 ) then
    ifault = 1
    return
  end if

  if ( n < 1 ) then
    ifault = 1
    return
  end if

  if ( konvge < 1 ) then
    ifault = 1
    return
  end if

  icount = 0
  numres = 0

  jcount = konvge
  dn = real ( n, kind = 8 )
  nn = n + 1
  dnn = real ( nn, kind = 8 )
  del = 1.0D+00
  rq = reqmin * dn
!
!  Initial or restarted loop.
!
  do

    do i = 1, n
      p(i,nn) = start(i)
    end do

    y(nn) = fn ( start )
    icount = icount + 1

    do j = 1, n
      x = start(j)
      start(j) = start(j) + step(j) * del
      do i = 1, n
        p(i,j) = start(i)
      end do
      y(j) = fn ( start )
      icount = icount + 1
      start(j) = x
    end do
!                    
!  The simplex construction is complete.
!                    
!  Find highest and lowest Y values.  YNEWLO = Y(IHI) indicates
!  the vertex of the simplex to be replaced.
!                
    ylo = y(1)
    ilo = 1

    do i = 2, nn
      if ( y(i) < ylo ) then
        ylo = y(i) 
        ilo = i
      end if
    end do
!
!  Inner loop.
!
    do

      if ( kcount <= icount ) then
        exit
      end if

      ynewlo = y(1)
      ihi = 1

      do i = 2, nn
        if ( ynewlo < y(i) ) then
          ynewlo = y(i)
          ihi = i
        end if
      end do
!
!  Calculate PBAR, the centroid of the simplex vertices
!  excepting the vertex with Y value YNEWLO.
!
      do i = 1, n
        z = 0.0D+00
        do j = 1, nn    
          z = z + p(i,j)
        end do
        z = z - p(i,ihi)   
        pbar(i) = z / dn   
      end do
!
!  Reflection through the centroid.
!
      do i = 1, n
        pstar(i) = pbar(i) + rcoeff * ( pbar(i) - p(i,ihi) )
      end do

      ystar = fn ( pstar )
      icount = icount + 1
!
!  Successful reflection, so extension.
!
      if ( ystar < ylo ) then

        do i = 1, n
          p2star(i) = pbar(i) + ecoeff * ( pstar(i) - pbar(i) )
        end do

        y2star = fn ( p2star )
        icount = icount + 1
!
!  Check extension.
!
        if ( ystar < y2star ) then

          do i = 1, n
            p(i,ihi) = pstar(i)
          end do

          y(ihi) = ystar
!
!  Retain extension or contraction.
!
        else

          do i = 1, n
            p(i,ihi) = p2star(i)
          end do

          y(ihi) = y2star

        end if
!
!  No extension.
!
      else

        l = 0
        do i = 1, nn
          if ( ystar < y(i) ) then
            l = l + 1
          end if
        end do

        if ( 1 < l ) then

          do i = 1, n
            p(i,ihi) = pstar(i)
          end do

          y(ihi) = ystar
!
!  Contraction on the Y(IHI) side of the centroid.
!
        else if ( l == 0 ) then

          do i = 1, n
            p2star(i) = pbar(i) + ccoeff * ( p(i,ihi) - pbar(i) )
          end do

          y2star = fn ( p2star )
          icount = icount + 1
!
!  Contract the whole simplex.
!
          if ( y(ihi) < y2star ) then

            do j = 1, nn
              do i = 1, n
                p(i,j) = ( p(i,j) + p(i,ilo) ) * 0.5D+00
                xmin(i) = p(i,j)
              end do
              y(j) = fn ( xmin )
              icount = icount + 1
            end do

            ylo = y(1)
            ilo = 1

            do i = 2, nn
              if ( y(i) < ylo ) then
                ylo = y(i) 
                ilo = i
              end if
            end do

            cycle
!
!  Retain contraction.
!
          else

            do i = 1, n
              p(i,ihi) = p2star(i)
            end do
            y(ihi) = y2star

          end if
!
!  Contraction on the reflection side of the centroid.
!
        else if ( l == 1 ) then
  
          do i = 1, n
            p2star(i) = pbar(i) + ccoeff * ( pstar(i) - pbar(i) )
          end do

          y2star = fn ( p2star )
          icount = icount + 1
!
!  Retain reflection?
!
          if ( y2star <= ystar ) then
 
            do i = 1, n
              p(i,ihi) = p2star(i)
            end do
            y(ihi) = y2star

          else

            do i = 1, n
              p(i,ihi) = pstar(i)
            end do
            y(ihi) = ystar  

          end if
 
        end if

      end if
!
!  Check if YLO improved.
!
      if ( y(ihi) < ylo ) then
        ylo = y(ihi)
        ilo = ihi
      end if

      jcount = jcount - 1

      if ( 0 < jcount ) then
        cycle
      end if
!
!  Check to see if minimum reached.
!
      if ( icount <= kcount ) then

        jcount = konvge

        z = 0.0D+00
        do i = 1, nn
          z = z + y(i)
        end do
        x = z / dnn

        z = 0.0D+00
        do i = 1, nn
          z = z + ( y(i) - x )**2
        end do

        if ( z <= rq ) then
          exit
        end if

      end if

    end do
!
!  Factorial tests to check that YNEWLO is a local minimum.
!
    do i = 1, n
      xmin(i) = p(i,ilo)
    end do

    ynewlo = y(ilo)

    if ( kcount < icount ) then
      ifault = 2
      exit
    end if

    ifault = 0

    do i = 1, n
      del = step(i) * eps
      xmin(i) = xmin(i) + del
      z = fn ( xmin )
      icount = icount + 1
      if ( z < ynewlo ) then
        ifault = 2
        exit
      end if
      xmin(i) = xmin(i) - del - del
      z = fn ( xmin )
      icount = icount + 1
      if ( z < ynewlo ) then
        ifault = 2
        exit
      end if
      xmin(i) = xmin(i) + del
    end do

    if ( ifault == 0 ) then
      exit
    end if
!
!  Restart the procedure.
!
    do i = 1, n
      start(i) = xmin(i)
    end do

    del = eps
    numres = numres + 1

  end do

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Modified:
!
!    06 August 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  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