subroutine fmin(fn,start,&
lambda,tol,conv_check,max_fun_calls,fun_calls,num_restart,ierr)
interface
function fn(x)
real(8),dimension(:) :: x
real(8) :: fn
end function fn
end interface
real(8) :: start(:)
real(8),optional :: lambda(size(start)) !--> step
real(8),optional :: tol !--> reqmin
integer,optional :: conv_check !--> konvge
integer,optional :: max_fun_calls !--> kcount
integer,optional :: fun_calls !--> icount
integer,optional :: num_restart !--> numres
integer,optional :: ierr !--> ifault
!
real(8) :: step(size(start))
real(8) :: reqmin
integer :: konvge
integer :: kcount
integer :: icount
integer :: numres
integer :: ifault
!
real(8) :: xmin(size(start))
integer :: n
real(8), parameter :: ccoeff = 0.5D+00
real(8) :: del
real(8), parameter :: ecoeff = 2.0D+00
real(8), parameter :: eps = 0.001D+00
integer :: i
integer :: ihi
integer :: ilo
integer :: j
integer :: jcount
integer :: l
real(8) :: p(size(start),size(start)+1)
real(8) :: p2star(size(start))
real(8) :: pbar(size(start))
real(8) :: pstar(size(start))
real(8),parameter :: rcoeff = 1.0D+00
real(8) :: rq
real(8) :: x
real(8) :: y(size(start)+1)
real(8) :: y2star
real(8) :: ylo
real(8) :: ynewlo
real(8) :: ystar
real(8) :: z
!
n = size(start)
!
step=1d0;if(present(lambda))step=lambda
reqmin=1d-8;if(present(tol))reqmin=tol
konvge=10;if(present(conv_check))konvge=conv_check
kcount=500;if(present(max_fun_calls))kcount=max_fun_calls
!
! 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
!
! Initialization.
!
icount = 0
numres = 0
jcount = konvge
del = 1.0D+00
rq = reqmin * real ( n, kind = 8 )
!
! Initial or restarted loop.
!
do
p(1:n,n+1) = start(1:n)
y(n+1) = fn ( start )
icount = icount + 1
!
! Define the initial simplex.
!
do j = 1, n
x = start(j)
start(j) = start(j) + step(j) * del
p(1:n,j) = start(1:n)
y(j) = fn ( start )
icount = icount + 1
start(j) = x
end do
!
! Find highest and lowest Y values. YNEWLO = Y(IHI) indicates
! the vertex of the simplex to be replaced.
!
ilo = minloc ( y(1:n+1), 1 )
ylo = y(ilo)
!
! Inner loop.
!
do while ( icount < kcount )
!
! YNEWLO is, of course, the HIGHEST value???
!
ihi = maxloc ( y(1:n+1), 1 )
ynewlo = y(ihi)
!
! Calculate PBAR, the centroid of the simplex vertices
! excepting the vertex with Y value YNEWLO.
!
do i = 1, n
pbar(i) = ( sum ( p(i,1:n+1) ) - p(i,ihi) ) / real ( n, kind = 8 )
end do
!
! Reflection through the centroid.
!
pstar(1:n) = pbar(1:n) + rcoeff * ( pbar(1:n) - p(1:n,ihi) )
ystar = fn ( pstar )
icount = icount + 1
!
! Successful reflection, so extension.
!
if ( ystar < ylo ) then
p2star(1:n) = pbar(1:n) + ecoeff * ( pstar(1:n) - pbar(1:n) )
y2star = fn ( p2star )
icount = icount + 1
!
! Retain extension or contraction.
!
if ( ystar < y2star ) then
p(1:n,ihi) = pstar(1:n)
y(ihi) = ystar
else
p(1:n,ihi) = p2star(1:n)
y(ihi) = y2star
end if
!
! No extension.
!
else
l = 0
do i = 1, n + 1
if ( ystar < y(i) ) then
l = l + 1
end if
end do
if ( 1 < l ) then
p(1:n,ihi) = pstar(1:n)
y(ihi) = ystar
!
! Contraction on the Y(IHI) side of the centroid.
!
else if ( l == 0 ) then
p2star(1:n) = pbar(1:n) + ccoeff * ( p(1:n,ihi) - pbar(1:n) )
y2star = fn ( p2star )
icount = icount + 1
!
! Contract the whole simplex.
!
if ( y(ihi) < y2star ) then
do j = 1, n + 1
p(1:n,j) = ( p(1:n,j) + p(1:n,ilo) ) * 0.5D+00
xmin(1:n) = p(1:n,j)
y(j) = fn ( xmin )
icount = icount + 1
end do
ilo = minloc ( y(1:n+1), 1 )
ylo = y(ilo)
cycle
!
! Retain contraction.
!
else
p(1:n,ihi) = p2star(1:n)
y(ihi) = y2star
end if
!
! Contraction on the reflection side of the centroid.
!
else if ( l == 1 ) then
p2star(1:n) = pbar(1:n) + ccoeff * ( pstar(1:n) - pbar(1:n) )
y2star = fn ( p2star )
icount = icount + 1
!
! Retain reflection?
!
if ( y2star <= ystar ) then
p(1:n,ihi) = p2star(1:n)
y(ihi) = y2star
else
p(1:n,ihi) = pstar(1:n)
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
x = sum ( y(1:n+1) ) / real ( n + 1, kind = 8 )
z = sum ( ( y(1:n+1) - x )**2 )
if ( z <= rq ) then
exit
end if
end if
end do
!
! Factorial tests to check that YNEWLO is a local minimum.
!
xmin(1:n) = p(1:n,ilo)
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.
!
start(1:n) = xmin(1:n)
del = eps
numres = numres + 1
end do
if(present(fun_calls))fun_calls=icount
if(present(num_restart))num_restart=numres
if(present(ierr))ierr=ifault
start = xmin
return
end subroutine fmin