199 lines
5.8 KiB
FortranFixed
199 lines
5.8 KiB
FortranFixed
subroutine cpRepeatCompassSearch(ndim,xbest,fbest,
|
|
& bmin,bmax,funkmin,f1dim,xtol)
|
|
implicit none
|
|
integer ndim
|
|
double precision xbest(1:ndim),fbest,
|
|
& bmin(1:ndim),bmax(1:ndim),xtol,f1dim
|
|
double precision fvalpre,dmax,xpre(1:ndim),ftol,direction(ndim)
|
|
parameter(ftol=1.0d-7)
|
|
integer i,n
|
|
logical resetran2
|
|
common /cpran2reset/resetran2
|
|
external funkmin,f1dim
|
|
!
|
|
n=0
|
|
resetran2=.true.
|
|
10 fvalpre=fbest
|
|
do i=1,ndim
|
|
xpre(i)=xbest(i)
|
|
enddo
|
|
call cpCompassSearch(ndim,xbest,fbest,
|
|
& bmin,bmax,funkmin,f1dim,xtol)
|
|
n=n+1
|
|
dmax=dabs(xbest(1)-xpre(1))
|
|
do i=2,ndim
|
|
if(dmax.lt.dabs(xbest(i)-xpre(i)))then
|
|
dmax=dabs(xbest(i)-xpre(i))
|
|
endif
|
|
enddo
|
|
if(dabs(fvalpre-fbest).gt.ftol.and.
|
|
& dmax.gt.xtol.and.n.lt.5000)then
|
|
do i=1,ndim
|
|
direction(i)=xbest(i)-xpre(i)
|
|
enddo
|
|
call linmin(xbest,bmin,bmax,direction,
|
|
& ndim,f1dim,fbest)
|
|
goto 10
|
|
endif
|
|
return
|
|
end subroutine cpRepeatCompassSearch
|
|
|
|
subroutine cpCompassSearch(ndim,xbest,fbest,
|
|
& bmin,bmax,funkmin,f1dim,xtol)
|
|
implicit none
|
|
|
|
! This subroutine minimizes the function funkmin using the compass search method. The ! maximum number of function evaluations is maxiter. Once mexiter is reached, all
|
|
! function evaluations are ranked and returned.
|
|
!
|
|
!------------------------------------- Inputs -----------------------------------------------------
|
|
! maxiter: the maximum number of function evaluations allowed
|
|
! xbest: the initial guess
|
|
! fbest: the cost function value at xinitial
|
|
! bmin: the lower bounds of the parameters to be optimized
|
|
! bmax: the upper bounds of the parameters to be optimized
|
|
! ndim: the number of parameters to optimize
|
|
! funkmin: the name of the function to minimize
|
|
|
|
!------------------------------------- Outputs ---------------------------------------------------
|
|
! xobs: points where the function is evaluated. Ranked from the best to worst with the
|
|
! first point being the best point.
|
|
! fvalue: the function values at xobs
|
|
! ierr: =0 convergence criterion not reached
|
|
! =1 convergence criterion reached (minimum found)
|
|
!
|
|
integer ndim
|
|
double precision xbest(1:ndim),fbest,f1dim,
|
|
& bmin(1:ndim),bmax(1:ndim),xtol,dx1,dx2
|
|
external funkmin,f1dim
|
|
!------------------------------- Locals -----------------------------------------------------------
|
|
double precision diftol,delta,
|
|
& xcompass(1:2*ndim,1:ndim),fcompass(1:2*ndim),
|
|
& xvec(1:ndim),xcent(1:ndim),fcent,dif,shrink,
|
|
& direction(ndim),dmax,fcent0,cpran2_reset
|
|
|
|
parameter(shrink=0.618d0,diftol=1.0d-7)
|
|
integer i,j,k
|
|
!
|
|
delta=0.618d0
|
|
do i=1,ndim
|
|
xcent(i)=xbest(i)
|
|
enddo
|
|
fcent=fbest
|
|
10 continue
|
|
do i=1,ndim
|
|
do j=1,ndim
|
|
xcompass(i,j)=xcent(j)
|
|
xcompass(ndim+i,j)=xcent(j)
|
|
enddo
|
|
xcompass(i,i)=xcent(i)+delta*(bmax(i)-xcent(i))
|
|
xcompass(ndim+i,i)=xcent(i)+delta*(bmin(i)-xcent(i))
|
|
enddo
|
|
|
|
do i=1,2*ndim
|
|
do j=1,ndim
|
|
xvec(j)=xcompass(i,j)
|
|
enddo
|
|
call funkmin(ndim,xvec,fcompass(i))
|
|
enddo
|
|
do i=1,ndim
|
|
xbest(i)=xcompass(1,i)
|
|
enddo
|
|
fbest=fcompass(1)
|
|
do i=2,2*ndim
|
|
if(fcompass(i).lt.fbest)then
|
|
fbest=fcompass(i)
|
|
do j=1,ndim
|
|
xbest(j)=xcompass(i,j)
|
|
enddo
|
|
endif
|
|
enddo
|
|
fcent0=fcent
|
|
do i=1,ndim
|
|
xvec(i)=xcent(i)
|
|
enddo
|
|
do i=1,ndim
|
|
dx1=xcompass(i,i)-xcent(i)
|
|
dx2=xcent(i)-xcompass(i+ndim,i)
|
|
direction(i)=0.0d0
|
|
if(dx1.ne.0.0d0)then
|
|
direction(i)=(fcompass(i)-fcent)/dx1
|
|
endif
|
|
if(dx2.ne.0.0d0)then
|
|
direction(i)=direction(i)+
|
|
& (fcent-fcompass(i+ndim))/dx2
|
|
endif
|
|
direction(i)=-0.5d0*direction(i)
|
|
if(direction(i).eq.0.0d0)direction(i)=
|
|
& cpran2_reset()-0.5d0
|
|
enddo
|
|
call cplinmin(xcent,bmin,bmax,direction,
|
|
& ndim,f1dim,fcent)
|
|
if(fcent.gt.fcent0)then
|
|
fcent=fcent0
|
|
do i=1,ndim
|
|
xcent(i)=xvec(i)
|
|
enddo
|
|
endif
|
|
dif=fcent-fbest
|
|
if(fbest.le.fcent)then
|
|
fcent=fbest
|
|
do i=1,ndim
|
|
xcent(i)=xbest(i)
|
|
enddo
|
|
endif
|
|
if(dif.ge.0.0d0)then
|
|
if(dif.gt.diftol)goto 10
|
|
if(delta.lt.diftol)goto 100
|
|
delta=delta*shrink
|
|
goto 10
|
|
else
|
|
!no progress
|
|
if(dabs(dif).gt.diftol)then
|
|
if(delta.lt.diftol)goto 100
|
|
delta=delta*shrink
|
|
goto 10
|
|
endif
|
|
dmax=dabs(xcompass(1,1)-xcompass(ndim+1,1))
|
|
do i=2,ndim
|
|
if(dmax.lt.dabs(xcompass(i,i)-
|
|
& xcompass(ndim+i,i)))then
|
|
dmax=dabs(xcompass(i,i)-
|
|
& xcompass(ndim+i,i))
|
|
endif
|
|
enddo
|
|
if(dmax.gt.xtol)then
|
|
if(delta.lt.diftol)goto 100
|
|
delta=delta*shrink
|
|
goto 10
|
|
else
|
|
goto 100
|
|
endif
|
|
endif
|
|
100 fbest=fcent
|
|
do i=1,ndim
|
|
xbest(i)=xcent(i)
|
|
dx1=xcompass(i,i)-xcent(i)
|
|
dx2=xcent(i)-xcompass(i+ndim,i)
|
|
direction(i)=0.0d0
|
|
if(dx1.ne.0.0d0)then
|
|
direction(i)=(fcompass(i)-fcent)/dx1
|
|
endif
|
|
if(dx2.ne.0.0d0)then
|
|
direction(i)=direction(i)+
|
|
& (fcent-fcompass(i+ndim))/dx2
|
|
endif
|
|
direction(i)=-0.5d0*direction(i)
|
|
if(direction(i).eq.0.0d0)direction(i)=
|
|
& cpran2_reset()-0.5d0
|
|
enddo
|
|
call cplinmin(xcent,bmin,bmax,direction,
|
|
& ndim,f1dim,fcent)
|
|
if(fcent.lt.fbest)then
|
|
fbest=fcent
|
|
do i=1,ndim
|
|
xbest(i)=xcent(i)
|
|
enddo
|
|
endif
|
|
return
|
|
end subroutine cpCompassSearch
|