122 lines
3.1 KiB
FortranFixed
122 lines
3.1 KiB
FortranFixed
SUBROUTINE amebsa(p,y,mp,np,ndim,pb,yb,ftol,funk,iter,temptr)
|
|
INTEGER iter,mp,ndim,np,NMAX
|
|
double precision ftol,temptr,yb,p(mp,np),pb(np),y(mp),funk
|
|
PARAMETER (NMAX=200)
|
|
EXTERNAL funk
|
|
CU USES amotsa,funk,ran1
|
|
INTEGER i,idum,ihi,ilo,inhi,j,m,n
|
|
double precision rtol,sum,swap,tt,yhi,ylo,ynhi,ysave,yt,ytry,
|
|
&psum(NMAX),amotsa,ran1
|
|
COMMON /ambsa/ tt,idum
|
|
tt=-temptr
|
|
1 do 12 n=1,ndim
|
|
sum=0.0d0
|
|
do 11 m=1,ndim+1
|
|
sum=sum+p(m,n)
|
|
11 continue
|
|
psum(n)=sum
|
|
12 continue
|
|
2 ilo=1
|
|
inhi=1
|
|
ihi=2
|
|
ylo=y(1)+tt*dlog(ran1(idum))
|
|
ynhi=ylo
|
|
yhi=y(2)+tt*dlog(ran1(idum))
|
|
if (ylo.gt.yhi) then
|
|
ihi=1
|
|
inhi=2
|
|
ilo=2
|
|
ynhi=yhi
|
|
yhi=ylo
|
|
ylo=ynhi
|
|
endif
|
|
do 13 i=3,ndim+1
|
|
yt=y(i)+tt*dlog(ran1(idum))
|
|
if(yt.le.ylo) then
|
|
ilo=i
|
|
ylo=yt
|
|
endif
|
|
if(yt.gt.yhi) then
|
|
inhi=ihi
|
|
ynhi=yhi
|
|
ihi=i
|
|
yhi=yt
|
|
else if(yt.gt.ynhi) then
|
|
inhi=i
|
|
ynhi=yt
|
|
endif
|
|
13 continue
|
|
rtol=2.0d0*dabs(yhi-ylo)/(dabs(yhi)+dabs(ylo))
|
|
if(rtol.lt.ftol.or.iter.lt.0) then
|
|
swap=y(1)
|
|
y(1)=y(ilo)
|
|
y(ilo)=swap
|
|
do 14 n=1,ndim
|
|
swap=p(1,n)
|
|
p(1,n)=p(ilo,n)
|
|
p(ilo,n)=swap
|
|
14 continue
|
|
return
|
|
endif
|
|
iter=iter-2
|
|
ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,-1.0d0)
|
|
if (ytry.le.ylo) then
|
|
ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,2.0d0)
|
|
else if (ytry.ge.ynhi) then
|
|
ysave=yhi
|
|
ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,0.5d0)
|
|
if (ytry.ge.ysave) then
|
|
do 16 i=1,ndim+1
|
|
if(i.ne.ilo)then
|
|
do 15 j=1,ndim
|
|
psum(j)=0.5d0*(p(i,j)+p(ilo,j))
|
|
p(i,j)=psum(j)
|
|
15 continue
|
|
y(i)=funk(psum)
|
|
endif
|
|
16 continue
|
|
iter=iter-ndim
|
|
goto 1
|
|
endif
|
|
else
|
|
iter=iter+1
|
|
endif
|
|
goto 2
|
|
END
|
|
|
|
double precision FUNCTION amotsa
|
|
&(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,fac)
|
|
INTEGER ihi,mp,ndim,np,NMAX
|
|
double precision amotsa,fac,yb,yhi,p(mp,np),pb(np),psum(np),
|
|
&y(mp),funk
|
|
PARAMETER (NMAX=200)
|
|
EXTERNAL funk
|
|
CU USES funk,ran1
|
|
INTEGER idum,j
|
|
double precision fac1,fac2,tt,yflu,ytry,ptry(NMAX),ran1
|
|
COMMON /ambsa/ tt,idum
|
|
fac1=(1.-fac)/ndim
|
|
fac2=fac1-fac
|
|
do 11 j=1,ndim
|
|
ptry(j)=psum(j)*fac1-p(ihi,j)*fac2
|
|
11 continue
|
|
ytry=funk(ptry)
|
|
if (ytry.le.yb) then
|
|
do 12 j=1,ndim
|
|
pb(j)=ptry(j)
|
|
12 continue
|
|
yb=ytry
|
|
endif
|
|
yflu=ytry-tt*log(ran1(idum))
|
|
if (yflu.lt.yhi) then
|
|
y(ihi)=ytry
|
|
yhi=yflu
|
|
do 13 j=1,ndim
|
|
psum(j)=psum(j)-p(ihi,j)+ptry(j)
|
|
p(ihi,j)=ptry(j)
|
|
13 continue
|
|
endif
|
|
amotsa=yflu
|
|
return
|
|
END
|