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,iter,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
integer, intent(out) :: iter
real(8), intent(in) :: eps
! local variables
integer i,j,il,iu
! 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,t1
! 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 : ",I8)') 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
iter=0
10 continue
iter=iter+1
if (iter >= maxit) return
! find the lowest and highest vertex
il=1
iu=1
do i=2,n+1
  if (f(i) < f(il)) il=i
  if (f(i) > f(iu)) iu=i
end do
! check for convergence
if ((f(iu)-f(il)) < eps) return
! compute the mean of the n lowest vertices
t1=1.d0/dble(n)
do i=1,n
  sm=0.d0
  do j=1,iu-1
    sm=sm+x(i,j)
  end do
  do j=iu+1,n+1
    sm=sm+x(i,j)
  end do
  xm(i)=t1*sm
end do
xr(:)=xm(:)+alpha*(xm(:)-x(:,iu))
fr=fmin_nm(id,rd,xr)
if (f(il) > fr) goto 30
if ((f(il) <= fr).and.(fr < f(iu))) then
! reflection
  x(:,iu)=xr(:)
  f(iu)=fr
  goto 10
else
  goto 40
end if
30 continue
xe(:)=xm(:)+gamma*(xr(:)-xm(:))
fe=fmin_nm(id,rd,xe)
if (fr > fe) then
! expansion
  x(:,iu)=xe(:)
  f(iu)=fe
else
! reflection
  x(:,iu)=xr(:)
  f(iu)=fr
end if
goto 10
40 continue
xc(:)=xm(:)+beta*(x(:,iu)-xm(:))
fc=fmin_nm(id,rd,xc)
if (fc < f(iu)) then
! contraction
  x(:,iu)=xc(:)
  f(iu)=fc
  goto 10
end if
! shrinkage
do j=1,il-1
  x(:,j)=x(:,il)+sigma*(x(:,j)-x(:,il))
  f(j)=fmin_nm(id,rd,x(1,j))
end do
do j=il+1,n+1
  x(:,j)=x(:,il)+sigma*(x(:,j)-x(:,il))
  f(j)=fmin_nm(id,rd,x(1,j))
end do
goto 10
end subroutine