404 lines
9.4 KiB
FortranFixed
404 lines
9.4 KiB
FortranFixed
SUBROUTINE powellann(p,xi,n,np,ftol,fret,pmin,pmax,
|
|
& funkmin,f1dim,ITMAX)
|
|
! fret must be given on entry
|
|
|
|
implicit none
|
|
INTEGER iter,n,np,NMAX,ITMAX
|
|
double precision fret,ftol,p(np),xi(np,np),TINY,
|
|
& pmin(np),pmax(np)
|
|
|
|
PARAMETER (NMAX=50,TINY=1.0d-25)
|
|
CU USES funkmin,annlinmin
|
|
INTEGER i,ibig,j
|
|
double precision del,fp,fptt,t,pt(NMAX),
|
|
& ptt(NMAX),xit(NMAX)
|
|
external funkmin,f1dim
|
|
|
|
do 11 j=1,n
|
|
pt(j)=p(j)
|
|
11 continue
|
|
iter=0
|
|
1 iter=iter+1
|
|
fp=fret
|
|
ibig=0
|
|
del=0.0d0
|
|
do 13 i=1,n
|
|
do 12 j=1,n
|
|
xit(j)=xi(j,i)
|
|
12 continue
|
|
fptt=fret
|
|
|
|
call annlinmin(p,pmin,pmax,xit,n,f1dim,fret)
|
|
if(fptt-fret.gt.del)then
|
|
del=fptt-fret
|
|
ibig=i
|
|
endif
|
|
13 continue
|
|
if(2.0d0*(fp-fret).le.ftol*(dabs(fp)+dabs(fret))+TINY)return
|
|
if(iter.eq.ITMAX)then
|
|
write(*,*)'powell exceeding maximum iterations'
|
|
return
|
|
endif
|
|
do 14 j=1,n
|
|
ptt(j)=2.0d0*p(j)-pt(j)
|
|
xit(j)=p(j)-pt(j)
|
|
pt(j)=p(j)
|
|
14 continue
|
|
call funkmin(n,ptt,fptt)
|
|
if(fptt.ge.fp)goto 1
|
|
t=2.0d0*(fp-2.0d0*fret+fptt)*(fp-fret-del)**2-
|
|
& del*(fp-fptt)**2
|
|
if(t.ge.0.0d0)goto 1
|
|
|
|
call annlinmin(p,pmin,pmax,xit,n,f1dim,fret)
|
|
|
|
do 15 j=1,n
|
|
xi(j,ibig)=xi(j,n)
|
|
xi(j,n)=xit(j)
|
|
15 continue
|
|
goto 1
|
|
END
|
|
C (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.
|
|
|
|
SUBROUTINE annlinmin(p,pmin,pmax,xi,n,f1dim,fret)
|
|
implicit none
|
|
INTEGER n,NMAX
|
|
double precision fret,p(n),xi(n),ftol,pmin(n),pmax(n)
|
|
PARAMETER (NMAX=1000,ftol=1.0d-6)
|
|
CU USES brent,f1dim,mnbrak
|
|
INTEGER j,ncom,k,ierr
|
|
double precision ax,bx,fa,fb,fx,xmin,xx,pcom(NMAX),
|
|
& xicom(NMAX),xxmin,xxmax,term,w(2),bph(2),
|
|
& paramnormsk,paramnormsb,ynorm,q(2),bend,terma,
|
|
& termb,termc,delta,reducer,root,f1dim,postannfunc
|
|
parameter(reducer=0.25d0)
|
|
|
|
COMMON /f1com/ pcom,xicom,ncom
|
|
|
|
integer maxm1dsamp,m1dsamp
|
|
parameter(maxm1dsamp=1000)
|
|
common /intannlinmin/m1dsamp
|
|
double precision y1dsamp(maxm1dsamp),params1d(maxm1dsamp)
|
|
common /dbleannlinmin/y1dsamp,params1d
|
|
|
|
EXTERNAL f1dim
|
|
|
|
ncom=n
|
|
delta=0.5d0
|
|
|
|
do j=1,n
|
|
pcom(j)=p(j)
|
|
xicom(j)=xi(j)
|
|
enddo
|
|
|
|
xxmax=1.0d+100
|
|
xxmin=-1.0d+100
|
|
do j=1,n
|
|
if(xicom(j).gt.1.0d-100)then
|
|
! if(xicom(j).gt.0.0d0)then
|
|
xx=(pmax(j)-pcom(j))/xicom(j)
|
|
ax=(pmin(j)-pcom(j))/xicom(j)
|
|
else
|
|
if(xicom(j).lt.(-1.0d-100))then
|
|
! if(xicom(j).lt.0.0d0)then
|
|
ax=(pmax(j)-pcom(j))/xicom(j)
|
|
xx=(pmin(j)-pcom(j))/xicom(j)
|
|
else
|
|
xx=1.0d+100
|
|
ax=-1.0d+100
|
|
endif
|
|
endif
|
|
if(xxmax.gt.xx)then
|
|
xxmax=xx
|
|
endif
|
|
if(xxmin.lt.ax)then
|
|
xxmin=ax
|
|
endif
|
|
enddo
|
|
ax=0.0d0
|
|
if(dabs(xxmax).gt.dabs(xxmin))then
|
|
xx=0.25d0*xxmax
|
|
else
|
|
xx=0.25d0*xxmin
|
|
endif
|
|
|
|
m1dsamp=0
|
|
call mnbrak(ax,xx,bx,fa,fx,fb,xxmin,xxmax,ierr,f1dim)
|
|
if(ierr.ne.0)then
|
|
! bracketing failed
|
|
fret=fx
|
|
do j=1,n
|
|
xi(j)=xx*xi(j)
|
|
p(j)=p(j)+xi(j)
|
|
enddo
|
|
return
|
|
endif
|
|
|
|
if(ax.gt.xx)then
|
|
xxmax=ax
|
|
xxmin=bx
|
|
else
|
|
xxmax=bx
|
|
xxmin=ax
|
|
endif
|
|
|
|
write(*,*)m1dsamp
|
|
|
|
if(m1dsamp.lt.6)then
|
|
term=xx+(xxmax-xx)*(delta+1.0d0)/2.0d0
|
|
fret=f1dim(term)
|
|
if(fret.lt.fx)then
|
|
xxmin=xx
|
|
xx=term
|
|
fx=fret
|
|
else
|
|
xxmax=term
|
|
endif
|
|
term=xx+(xxmax-xx)*(1.0d0-delta)/2.0d0
|
|
fret=f1dim(term)
|
|
if(fret.lt.fx)then
|
|
xxmin=xx
|
|
xx=term
|
|
fx=fret
|
|
else
|
|
xxmax=term
|
|
endif
|
|
term=xx-(xx-xxmin)*(1.0d0+delta)/2.0d0
|
|
fret=f1dim(term)
|
|
if(fret.lt.fx)then
|
|
xxmax=xx
|
|
xx=term
|
|
fx=fret
|
|
else
|
|
xxmin=term
|
|
endif
|
|
term=xx-(xx-xxmin)*(1.0d0-delta)/2.0d0
|
|
fret=f1dim(term)
|
|
if(fret.lt.fx)then
|
|
xxmax=xx
|
|
xx=term
|
|
fx=fret
|
|
else
|
|
xxmin=term
|
|
endif
|
|
endif
|
|
|
|
w(1)=0.75334d0
|
|
w(2)=0.13425d0
|
|
bph(1)=0.01d0
|
|
bph(2)=-0.07d0
|
|
q(1)=1.2d0
|
|
q(2)=-2.0d0
|
|
|
|
100 call annfitting1d(m1dsamp,params1d,y1dsamp,paramnormsk,
|
|
& paramnormsb,ynorm,w,bph,q,bend)
|
|
|
|
|
|
term=q(1)*w(1)
|
|
fret=q(2)*w(2)
|
|
terma=term*w(2)*w(2)+fret*w(1)*w(1)
|
|
termb=2.0d0*(term*w(2)*bph(2)+fret*w(1)*bph(1))
|
|
termc=term*(1.0d0+bph(2)*bph(2))+
|
|
& fret*(1.0d0+bph(1)*bph(1))
|
|
|
|
if(terma.eq.0.0d0)then
|
|
if(termb.eq.0.0d0)goto 200
|
|
root=-termc/termb
|
|
goto 110
|
|
endif
|
|
fret=termb*termb-4.0d0*terma*termc
|
|
if(fret.lt.0.0d0)goto 200
|
|
term=-0.5d0*(termb+dsign(1.0d0,termb)*dsqrt(fret))
|
|
root=termc/term
|
|
root=(root-paramnormsb)/paramnormsk
|
|
if(root.lt.xxmax.and.root.gt.xxmin)then
|
|
if(dabs(root-xx).lt.ftol)goto 1000
|
|
fret=f1dim(root)
|
|
if(fret.lt.fx)then
|
|
if(root.lt.xx)then
|
|
xxmax=xx
|
|
else
|
|
xxmin=xx
|
|
endif
|
|
xx=root
|
|
fx=fret
|
|
goto 100
|
|
else
|
|
if(root.lt.xx)then
|
|
xxmin=root
|
|
else
|
|
xxmax=root
|
|
endif
|
|
endif
|
|
endif
|
|
root=term/terma
|
|
110 root=(root-paramnormsb)/paramnormsk
|
|
if(root.lt.xxmax.and.root.gt.xxmin)then
|
|
if(dabs(root-xx).lt.ftol)goto 1000
|
|
fret=f1dim(root)
|
|
if(fret.lt.fx)then
|
|
if(root.lt.xx)then
|
|
xxmax=xx
|
|
else
|
|
xxmin=xx
|
|
endif
|
|
xx=root
|
|
fx=fret
|
|
goto 100
|
|
else
|
|
if(root.lt.xx)then
|
|
xxmin=root
|
|
else
|
|
xxmax=root
|
|
endif
|
|
endif
|
|
endif
|
|
|
|
200 term=xx+delta*(xxmax-xx)
|
|
|
|
write(*,*)'ann failed'
|
|
|
|
fret=f1dim(term)
|
|
if(fret.lt.fx)then
|
|
xxmin=xx
|
|
fx=fret
|
|
xx=term
|
|
goto 100
|
|
endif
|
|
xxmax=term
|
|
term=xx-delta*(xx-xxmin)
|
|
fret=f1dim(term)
|
|
if(fret.lt.fx)then
|
|
xxmax=xx
|
|
fx=fret
|
|
xx=term
|
|
goto 100
|
|
endif
|
|
xxmin=term
|
|
delta=delta*reducer
|
|
if(delta.gt.ftol)goto 200
|
|
|
|
1000 fret=fx
|
|
do j=1,n
|
|
xi(j)=xx*xi(j)
|
|
p(j)=p(j)+xi(j)
|
|
enddo
|
|
|
|
return
|
|
END
|
|
C (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.
|
|
|
|
SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,xxmin,xxmax,
|
|
& ierr,func)
|
|
implicit none
|
|
double precision ax,bx,cx,fa,fb,fc,
|
|
& func,GOLD,GLIMIT,TINY
|
|
EXTERNAL func
|
|
PARAMETER(GOLD=1.618034d0,GLIMIT=100.0d0,TINY=1.0d-20)
|
|
double precision dum,fu,q,r,u,ulim,xxmin,xxmax
|
|
integer ierr
|
|
ierr=0
|
|
fa=func(ax)
|
|
fb=func(bx)
|
|
if(fb.gt.fa)then
|
|
dum=ax
|
|
ax=bx
|
|
bx=dum
|
|
dum=fb
|
|
fb=fa
|
|
fa=dum
|
|
endif
|
|
if(fa.eq.fb)then
|
|
cx=(bx+ax)/2.0d0
|
|
fc=func(cx)
|
|
if(fc.le.fa)return
|
|
endif
|
|
cx=bx+GOLD*(bx-ax)
|
|
if(cx.le.xxmin)then
|
|
cx=0.5d0*(dmin1(ax,bx)+xxmin)
|
|
endif
|
|
if(cx.ge.xxmax)then
|
|
cx=0.5d0*(dmax1(ax,bx)+xxmax)
|
|
endif
|
|
fc=func(cx)
|
|
1 if(fb.ge.fc)then
|
|
r=(bx-ax)*(fb-fc)
|
|
q=(bx-cx)*(fb-fa)
|
|
u=bx-((bx-cx)*q-(bx-ax)*r)/
|
|
& (2.0d0*dsign(dmax1(dabs(q-r),TINY),q-r))
|
|
ulim=bx+GLIMIT*(cx-bx)
|
|
if(ulim.ge.xxmax)then
|
|
ulim=xxmax-tiny
|
|
endif
|
|
if(ulim.le.xxmin)then
|
|
ulim=xxmin+tiny
|
|
endif
|
|
if((bx-u)*(u-cx).gt.0.0d0)then
|
|
fu=func(u)
|
|
if(fu.lt.fc)then
|
|
ax=bx
|
|
fa=fb
|
|
bx=u
|
|
fb=fu
|
|
return
|
|
elseif(fu.gt.fb)then
|
|
cx=u
|
|
fc=fu
|
|
return
|
|
endif
|
|
u=cx+GOLD*(cx-bx)
|
|
if(u.gt.xxmax)then
|
|
u=cx+0.5d0*(xxmax-cx)
|
|
endif
|
|
if(u.lt.xxmin)then
|
|
u=cx+0.5d0*(xxmin-cx)
|
|
endif
|
|
fu=func(u)
|
|
elseif((cx-u)*(u-ulim).gt.0.0d0)then
|
|
fu=func(u)
|
|
if(fu.lt.fc)then
|
|
bx=cx
|
|
cx=u
|
|
u=cx+GOLD*(cx-bx)
|
|
if(u.gt.xxmax)then
|
|
u=cx+0.5d0*(xxmax-cx)
|
|
endif
|
|
if(u.lt.xxmin)then
|
|
u=cx+0.5d0*(xxmin-cx)
|
|
endif
|
|
fb=fc
|
|
fc=fu
|
|
fu=func(u)
|
|
endif
|
|
else if((u-ulim)*(ulim-cx).ge.0.0d0)then
|
|
u=ulim
|
|
fu=func(u)
|
|
else
|
|
u=cx+GOLD*(cx-bx)
|
|
if(u.gt.xxmax)then
|
|
u=cx+0.5d0*(xxmax-cx)
|
|
endif
|
|
if(u.lt.xxmin)then
|
|
u=cx+0.5d0*(xxmin-cx)
|
|
endif
|
|
fu=func(u)
|
|
endif
|
|
ax=bx
|
|
bx=cx
|
|
cx=u
|
|
fa=fb
|
|
fb=fc
|
|
fc=fu
|
|
r=dmin1(dabs(ax-bx),dabs(ax-cx))
|
|
r=dmin1(r,dabs(bx-cx))
|
|
if(r.lt.tiny)then
|
|
! bracketing failed
|
|
ierr=1
|
|
return
|
|
endif
|
|
goto 1
|
|
endif
|
|
return
|
|
END
|
|
C (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.
|