! ----------------------------------------------------------- ! ! CONJUGATE GRADIENT MINIMIZATION MODULE ADOPTED FROM ! ! NUMERICAL RECIPES. REWRITTEN IN F90. ! ! ----------------------------------------------------------- ! MODULE CG_MINIMIZER ! ======================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER, PRIVATE :: NOP, NCOM DOUBLE PRECISION, PRIVATE, ALLOCATABLE, DIMENSION(:) :: PCOM,XICOM DOUBLE PRECISION, PRIVATE, ALLOCATABLE, DIMENSION(:) :: G,H,XI CONTAINS SUBROUTINE CG_INITIALIZE (NUMBER_OF_OPERAND) ! ---------------------------------------------------------------------- INTEGER NUMBER_OF_OPERAND NOP = NUMBER_OF_OPERAND ALLOCATE (PCOM(NOP)) ALLOCATE (XICOM(NOP)) ALLOCATE (G(NOP)) ALLOCATE (H(NOP)) ALLOCATE (XI(NOP)) END SUBROUTINE CG_INITIALIZE ! ---------------------------------------------------------------------- SUBROUTINE CG_FINALIZE ! ---------------------------------------------------------------------- DEALLOCATE (PCOM) DEALLOCATE (XICOM) DEALLOCATE (G) DEALLOCATE (H) DEALLOCATE (XI) END SUBROUTINE CG_FINALIZE ! ---------------------------------------------------------------------- SUBROUTINE MNBRAK(AX,BX,CX,FA,FB,FC,FUNC) ! ---------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) EXTERNAL func double precision ax,bx,cx,fa,fb,fc,func double precision, parameter :: & GOLD=1.618034,GLIMIT=1000.d0,TINY=1.d-20 double precision dum,fu,q,r,u,ulim fa=func(pcom+ax*xicom) fb=func(pcom+bx*xicom) if(fb.gt.fa)then dum=ax ax=bx bx=dum dum=fb fb=fa fa=dum endif cx=bx+GOLD*(bx-ax) fc=func(pcom+cx*xicom) 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.d0*sign(max(abs(q-r),TINY),q-r)) ulim=bx+GLIMIT*(cx-bx) if((bx-u)*(u-cx).gt.0.d0)then fu=func(pcom+u*xicom) if(fu.lt.fc)then ax=bx fa=fb bx=u fb=fu return else if(fu.gt.fb)then cx=u fc=fu return endif u=cx+GOLD*(cx-bx) fu=func(pcom+u*xicom) else if((cx-u)*(u-ulim).gt.0.d0)then fu=func(pcom+u*xicom) if(fu.lt.fc)then bx=cx cx=u u=cx+GOLD*(cx-bx) fb=fc fc=fu fu=func(pcom+u*xicom) endif else if((u-ulim)*(ulim-cx).ge.0.d0)then u=ulim fu=func(pcom+u*xicom) else u=cx+GOLD*(cx-bx) fu=func(pcom+u*xicom) endif ax=bx bx=cx cx=u fa=fb fb=fc fc=fu goto 1 endif END SUBROUTINE MNBRAK ! ---------------------------------------------------------------------- FUNCTION BRENT(AX,BX,CX,FUNC,TOL,XMIN) ! ---------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) integer, parameter :: itmax=1000 double precision, parameter :: CGOLD=.3819660, ZEPS=1.0d-20 external func double precision ax,bx,cx,tol,xmin,func, BRENT integer iter double precision a,b,d,e,etemp,fu,fv,fw,fx, & p,q,r,tol1,tol2,u,v,w,x,xm a=min(ax,cx) b=max(ax,cx) v=bx w=v x=v e=0.d0 fx=func(pcom+x*xicom) fv=fx fw=fx do 11 iter=1,ITMAX xm=0.5*(a+b) tol1=tol*abs(x)+ZEPS tol2=2.d0*tol1 if(abs(x-xm).le.(tol2-.5*(b-a))) goto 3 if(abs(e).gt.tol1) then r=(x-w)*(fx-fv) q=(x-v)*(fx-fw) p=(x-v)*q-(x-w)*r q=2.*(q-r) if(q.gt.0.d0) p=-p q=abs(q) etemp=e e=d if(abs(p).ge.abs(.5*q*etemp).or.p.le.q*(a-x).or.p.ge.q*(b-x)) goto 1 d=p/q u=x+d if(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x) goto 2 endif 1 if(x.ge.xm) then e=a-x else e=b-x endif d=CGOLD*e 2 if(abs(d).ge.tol1) then u=x+d else u=x+sign(tol1,d) endif fu=func(pcom+u*xicom) if(fu.le.fx) then if(u.ge.x) then a=x else b=x endif v=w fv=fw w=x fw=fx x=u fx=fu else if(u.lt.x) then a=u else b=u endif if(fu.le.fw .or. w.eq.x) then v=w fv=fw w=u fw=fu else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then v=u fv=fu endif endif 11 continue stop 'brent exceed maximum iterations' 3 xmin=x brent=fx END FUNCTION BRENT ! ---------------------------------------------------------------------- SUBROUTINE LINMIN (P,myxi,N,FRET,VALUE) ! ---------------------------------------------------------------------- integer n external value double precision fret,p(:),myxi(:),value ! tolerance passed to brent double precision, parameter :: tol=2.0d-4 double precision ax,bx,fa,fb,fx,xmin,xx ncom=n pcom=p xicom=myxi ax=0.d0 xx=1.d0 call mnbrak(ax,xx,bx,fa,fx,fb,VALUE) fret=brent(ax,xx,bx,VALUE,TOL,xmin) myxi=xmin*myxi p=p+myxi END SUBROUTINE LINMIN ! ----------------------------------------------------------------------- SUBROUTINE CG_MINIMIZE (P,N,FTOL,ITER,FRET,VALUE,FORCE) ! ----------------------------------------------------------------------- integer iter,n external value, force double precision fret,ftol,p(:),value,fp,dgg,gg integer, parameter :: ITMAX=10000 double precision, parameter :: EPS=1.d-20 iter=0 fp=value(p) call force(p,xi) if(dot_product(xi,xi).eq.0.d0) return g=-xi h=g xi=h do 14 iter=1,ITMAX call linmin(p,xi,n,fret,value) if(2.*abs(fret-fp).le.ftol*(abs(fret)+abs(fp)+EPS)) return fp=value(p) call force(p,xi) gg=dot_product(g,g) dgg=dot_product(xi+g,xi) if(gg.eq.0.d0) return gam=dgg/gg g=-xi h=g+gam*h xi=h 14 continue stop 'frprmn maximum iterations exceeded' END SUBROUTINE CG_MINIMIZE ! ----------------------------------------------------------------------- END MODULE CG_MINIMIZER ! =======================================================================