minf_nm.f90

Go to the documentation of this file.

! Copyright (C) 2007 J. K. Dewhurst and D. W. H. Rankin.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine minf_nm(id,rd,n,x,maxit,eps)
implicit none
! arguments
integer, intent(in) :: id(*)
real(8), intent(in) :: rd(*)
integer, intent(in) :: n
real(8), intent(inout) :: x(n,n+1)
integer, intent(in) :: maxit
real(8), intent(in) :: eps
! local variables
integer it,il,iu,i,j
! Nelder-Mead parmeters
real(8), parameter :: alpha=1.d0,gamma=2.d0
real(8), parameter :: beta=0.5d0,sigma=0.5d0
real(8) fr,fe,fc,sm
! automatic arrays
real(8) f(n+1),xm(n),xr(n),xe(n),xc(n)
! external functions
real(8), external :: fmin_nm
if (n < 1) then
  write(*,*)
  write(*,'("Error(minf_nm): n < 1 : ",I0)') n
  write(*,*)
  stop
end if
! evaluate the function at each vertex
do i=1,n+1
  f(i)=fmin_nm(id,rd,x(:,i))
end do
! begin iteration loop
do it=1,maxit
! find the lowest and highest vertex
  il=minloc(f(1:n+1),1)
  iu=maxloc(f(1:n+1),1)
! check for convergence
  if ((f(iu)-f(il)) < eps) return
! compute the mean of the n lowest vertices
  do i=1,n
    sm=sum(x(i,1:iu-1))+sum(x(i,iu+1:n+1))
    xm(i)=sm/n
  end do
  xr(1:n)=xm(1:n)+alpha*(xm(1:n)-x(1:n,iu))
  fr=fmin_nm(id,rd,xr)
  if (f(il) < fr) then
    if (fr < f(iu)) then
! reflection
      x(1:n,iu)=xr(1:n)
      f(iu)=fr
    else
      xc(1:n)=xm(1:n)+beta*(x(1:n,iu)-xm(1:n))
      fc=fmin_nm(id,rd,xc)
      if (fc < f(iu)) then
! contraction
        x(1:n,iu)=xc(1:n)
        f(iu)=fc
      else
! shrinkage
        do j=1,n+1
          if (j == il) cycle
          x(1:n,j)=x(1:n,il)+sigma*(x(1:n,j)-x(1:n,il))
          f(j)=fmin_nm(id,rd,x(1,j))
        end do
      end if
    end if
  else
    xe(1:n)=xm(1:n)+gamma*(xr(1:n)-xm(1:n))
    fe=fmin_nm(id,rd,xe)
    if (fr > fe) then
! expansion
      x(1:n,iu)=xe(1:n)
      f(iu)=fe
    else
! reflection
      x(1:n,iu)=xr(1:n)
      f(iu)=fr
    end if
  end if
end do
end subroutine