cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Nov 7 2000 c Dez 8 2000 cc Versions c Feb. 2001 c May 2001 c July 2001 cccccccccccccccccccccccccccccccccccccccccccccc c RGF modular version c cccccccccccccccccccccccccccccccccccccccccccccc c c Revision 2.2, Aug 31 for TASC-adaption, c note: utils.f also are changed c c Michael Hirsch / Wolfgang Quapp c Uni. Leipzig cc Block Data files integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene common /iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene data isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene c /20,21,25,31,32,33,34,35,36/ end program rgf real*8 crit_stop,crit_pc,pstepl,some_real,eold integer nn,irgfmain,istat,itmin,itmax,irgf_para,narg,iargc integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene, c job_id,job,igstep character*40 z_FoP,z_FoG,z_FoH,z_FoB,z_FoR,z_store,z_para, c z_msg,z_FoE logical ostore,ostat,o_flat,o_free,o_msg,o_no_up,o_const_sgn, c o_even,o_recorr,o_nn,o_corr,o_req_ex,o_SL dimension ostat(8),irgf_para(5),some_real(8) common /iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene common /iofiles/z_store,z_para,z_msg,z_FoR,z_FoP, c z_FoG,z_FoH,z_FoB,z_FoE common /para/crit_stop,crit_pc,pstepl,itmin,itmax,ostat common /old/job_id,irgf_para,o_corr,o_req_ex,o_SL,igstep,eold common /switch/o_flat,o_free,o_msg,o_no_up,o_const_sgn, c o_even,o_recorr,o_nn namelist /init/job_id,nn namelist /param/job_id,nn,z_FoP,z_FoG,z_FoH,z_FoB,z_FoR, c o_flat,o_free,o_msg,o_no_up,o_const_sgn,o_even,o_recorr,o_nn, c crit_stop,crit_pc,pstepl,itmin,itmax,z_FoE namelist /STORE/job,irgf_para,pstepl,crit_pc,crit_stop, c o_corr,o_req_ex,o_SL,igstep,eold,some_real data z_store/'store.rgf'/, z_para/'param.rgf'/, c z_msg/'message.rgf'/ c command line arguments c switches will be stored in OSTAT(3)...OSTAT(8) narg=iargc() if (narg.gt.0) call handleargs(narg,ostat) inquire(file=z_para,exist=ostat(1),opened=ostat(2)) if (.not.ostat(1)) call exit(1) if (ostat(2)) close(ipar) open(ipar,file=z_para,status='old') read(ipar,param,iostat=istat) if (istat.ne.0.and.istat.ne.119) then write(*,*)'wrong format of parameter file' call help(6,1) call exit(1) endif close(ipar) inquire(file=z_store,exist=ostat(1)) if (ostat(1)) then open(isto,file=z_store,status='old',iostat=istat) read(isto,nml=STORE,iostat=istat,err=1111) else job=job_id+1 endif if (job.eq.job_id) then ostore=.true. else ostore=.false. if (o_msg) call recreate(imsg,z_msg) if (o_msg) call recreate(ires,z_FoR) job=job_id endif if (o_msg) then open(imsg,file=z_msg,status='unknown', c access='append') if (ostore) then write(imsg,10) else write(imsg,20) endif write(imsg,*)'job_id Id:',job_id endif call closys istat=irgfmain(nn,ostore) call exit (istat) 1111 continue write(*,*) 'NameList Lesefehler',istat call exit cmd return 10 format(1x,'- use old "store.rgf" file. ') 20 format(1x,'- new job. create new "store.rgf"') end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c C main routine C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC integer function irgfmain(nn,ostore) logical ostore,oerror,omat_inv,oup,o_corr,o_bif,o_turn,o_req_ex logical o_flat,o_free,o_msg,o_no_up,o_const_sgn, c o_even,o_recorr,o_nn,ostat,o_SL integer i,j,nn,itmin,itmax,irgf_para,iochk, c istep,ipstep,isig,i_dfp,igstep integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene,job real*8 crit_stop,crit_pc,dsqrt,hessian,cold,fstep, c rsel,bmat,dk,bpinv,gmat,ginv,unit,rgeq,sw,rgrad,point,val,adj, c gradient,gradcont,tang,tcov,t1,t1old,told,pstepl,gnorm,dnst, c vec_ddot,drsel,dstep,rgnorm,coef,detk,det,tgt,tg,newton,dmnorm, c eival,eivec,diag,speinv,dxy,dx,dg,hold,xold,gold,tgf,proj, c energy,eold c real*8 some_real,ray,snorm,dray,dsn,bn,bnorm character*40 z_store,z_para,z_msg,z_FoR,z_FoP, c z_FoG,z_FoH,z_FoB,z_FoE character*20 zf_prop character*80 zval dimension hessian(nn,nn),gradient(nn),fstep(nn) dimension rsel(nn),bmat(nn,nn+6),dk(nn,nn+1), c bpinv(nn+6,nn),rgeq(nn-1,nn),dstep(nn),proj(nn-1,nn) dimension sw(nn),rgrad(nn-1),point(nn), c gradcont(nn),eival(nn),eivec(nn,nn),isig(3),val(4) dimension tcov(nn),tang(nn),t1(nn),t1old(nn),told(nn), c diag(nn,nn),speinv(nn,nn),dxy(nn,nn),cold(nn),irgf_para(5), c gmat(nn,nn),ginv(nn,nn),unit(nn,nn),adj(nn,nn) dimension hold(nn,nn),gold(nn),xold(nn),dx(nn),dg(nn),some_real(8) dimension ostat(8) common /iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene common /iofiles/z_store,z_para,z_msg,z_FoR,z_FoP, c z_FoG,z_FoH,z_FoB,z_FoE common /para/crit_stop,crit_pc,pstepl,itmin,itmax,ostat common /old/job,irgf_para,o_corr,o_req_ex,o_SL,igstep,eold common /switch/o_flat,o_free,o_msg,o_no_up,o_const_sgn, c o_even,o_recorr,o_nn namelist /store/job,irgf_para,pstepl,crit_pc,crit_stop, c o_corr,o_req_ex,some_real data zval(1:40)/'determinant of K (tang-2,told-2) '/ data zval(41:80)/'(tang-2,firststep) (tang-1,told-1) '/ c c Initialisierung und Einlesen einiger Paramter o_bif=.false. o_turn=.false. call vec_dinit(some_real,8,0.0d0) open(imsg,file=z_msg,status='unknown',access='append') open(ipar,file=z_para,status='old',iostat=iochk) c call openext call mat_diag(unit,nn,1.0d0) if (ostore) then open(isto,file=z_store,status='old') read(isto,*)fstep read(isto,*)proj read(isto,*)told read(isto,*)t1old read(isto,*)cold read(isto,*)xold read(isto,*)gold read(isto,*)hold read(isto,*)eold read(isto,nml=store) if (ostat(5)) call reread(z_para) close(isto) c irgf_para: c 1 - istep c 2 - ipstep c 3 - itmin c 4 - itmax c 5 - signature of the hessian c 6 - exact hessian next point (1 yes, 0 no) c 7 - corrector this point (1 yes, 0 predictor) istep=irgf_para(1) ipstep=irgf_para(2) itmin=irgf_para(3) itmax=irgf_para(4) else call vec_iinit(irgf_para,5,0) ipstep=0 istep=0 if (.not.o_free) then rewind(ipar) read (ipar,*)rsel oerror=.true. do i=1,nn if (rsel(i).ne.0.0d0) oerror=.false. enddo if (oerror) call err(-20,'rsel',3) call ortcompl(rsel,proj,unit,nn) endif endif call openext read(ipnt,*,iostat=iochk)point read(igrd,*,iostat=iochk)gradient oup=(istep.gt.0).and.(.not.o_no_up).and.(.not.o_req_ex) c force exact vs. updated HESSIAN by command line args oup=(oup.or.ostat(4)).and.(.not.ostat(3)) if (.not.oup) then read(ihes,*,iostat=iochk)hessian call message('(1x,''=== exact Hessian ==='')',28) endif if (.not.o_flat) c read(iwil,*,iostat=iochk)bmat read(iene,*,iostat=iochk)energy cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c if (o_flat) then call vec_dcopy(unit,gmat,nn*nn) call vec_dcopy(unit,ginv,nn*nn) else call getmet(bmat,bpinv,ginv,gmat,nn) endif if ((istep.eq.0).and.(o_free)) then call matmult(ginv,nn,gradient,1,rsel,nn) call ortcompl(rsel,proj,unit,nn) endif if (oup) then call vec_dmult_add(gradient,gold,nn,-1.0d0,dg) call vec_dmult_add(point,xold,nn,-1.0d0,dx) call err(i_dfp(hold,dg,dx,nn,hessian),'DFP',3) call message('(1x,''=== updated Hessian ==='')',30) endif call vec_dcopy(gradient,gold,nn) call vec_dcopy(point,xold,nn) call vec_dcopy(hessian,hold,nn*nn) call toang(point,sw,nn) call nka(ginv,hessian,eival,eivec,nn) call sig(eival,isig,nn) j=100*isig(1)+isig(2) if ((j.ne.irgf_para(5)).and.o_const_sgn) then o_req_ex=.true. call message c ('(1x,''--- exact hessian in the next step ---'')',45) else o_req_ex=.false. endif c signature of hessian irgf_para(5)=j 776 format(1x,'signature of hessian: ',i2,'+ ',i2,'- ',i2,'=') if (o_msg) write(imsg,776)(isig(i),i=1,3) call vec_dinit(diag,nn*nn,0.0d0) do i=1,nn diag(i,i)=eival(i) enddo if (omat_inv(eivec,nn,speinv)) then call matmult(diag,nn,gmat,nn,dxy,nn) endif c if (.not.ostore) call ortcompl(rsel,proj,gmat,nn) dnst=newton(hessian,gradient,dstep,gmat,nn) call vec_dinit(gradcont,nn,0.0d0) call blas_dgemm $ ('n','n',nn,1,nn,1.0d0,ginv, $ nn,gradient,nn,1.0d0,gradcont,nn) gnorm=dsqrt(vec_ddot(gradient,gradcont,nn)) if (.not.ostore) then write(imsg,*)' - calculation of first step' drsel=dmnorm(rsel,gmat,nn) if (o_even.and.o_free) then call vec_dmult_constant(rsel,nn,-(pstepl/drsel),fstep) else call vec_dmult_constant(rsel,nn,pstepl/drsel,fstep) endif endif istep=istep+1 coef=1.0d-5 if (istep.eq.1) then call vec_dcopy(fstep,dstep,nn) call mat_adj(hessian,nn,adj) o_corr=.false. call matmult(fstep,1,gmat,nn,tcov,nn) call vec_dcopy(fstep,t1,nn) call matmult(fstep,1,gmat,nn,t1old,nn) call matmult(fstep,1,gmat,nn,told,nn) else call tri_mat(proj,nn-1,nn,ginv,nn,hessian,nn,rgeq) call vec_dinit(rgrad,nn-1,0.0d0) call matmult(proj,nn-1,gradcont,1,rgrad,nn) call vec_dinit(sw,nn,0.0d0) do i=1,nn-1 sw(i+1)=rgrad(i) enddo cccccccccccccccccccccc rgnorm=dmnorm(sw,gmat,nn) call gettang(rgeq,gmat,tang,nn) call vec_dinit(sw,nn,0.0d0) cccccccccccccccccccccc call matmult(rgeq,nn-1,tang,1,sw,nn) call vec_dinit(tcov,nn,0.0d0) call blas_dgemm * ('n','n',nn,1,nn,1.0d0,gmat, * nn,tang,nn,1.0d0,tcov,nn) call buildk(dk,rgeq,tcov,rgrad,nn) detk=det(dk,nn) if (detk.lt.0.0d0) then call vec_dmult_constant(tang,nn,-1.0d0,t1) else call vec_dcopy(tang,t1,nn) endif tgt=vec_ddot(tang,told,nn) if (tgt.lt.0.0d0) then call vec_dmult_constant(tang,nn,-1.0d0,tang) call vec_dmult_constant(tcov,nn,-1.0d0,tcov) endif call buildk(dk,rgeq,tcov,rgrad,nn) call mat_adj(hessian,nn,adj) if (rgnorm.gt.crit_pc) then c Corrector call linsolve(dk,dstep,nn) if (o_recorr) c call recorrect(dstep,cold,gmat,o_corr,nn) if (o_corr) call message('(1x,''recorrect'')',16) o_corr=.true. else c Predictor ipstep=ipstep+1 call vec_dmult_constant(tang,nn,pstepl,dstep) if (ipstep.gt.1) then tgf=vec_ddot(fstep,tcov,nn) tg=vec_ddot(t1,t1old,nn) if (tg.lt.0.0d0) o_bif=.true. if (tgf.lt.0.0d0) then o_turn=.true. call vec_dmult_constant(fstep,nn,-1.0d0,fstep) endif endif call vec_dcopy(tcov,told,nn) call vec_dinit(t1old,nn,0.0d0) call blas_dgemm('n','n',nn,1,nn,1.0d0,gmat, * nn,t1,nn,1.0d0,t1old,nn) o_corr=.false. dsn=snorm(adj,gradient,ginv,gmat,nn) dray=ray(adj,gradient,ginv,nn) call func_prop * (some_real(1),some_real(2),dsn,iochk,zf_prop) write(imsg,1140)dsn,zf_prop call func_prop * (some_real(3),some_real(4),dray,iochk,zf_prop) write(imsg,1150)dray,zf_prop some_real(1)=some_real(2) some_real(2)=dsn some_real(3)=some_real(4) some_real(4)=dray if (dsn.gt.some_real(5)) some_real(5)=dsn endif endif val(1)=detk val(2)=tgt val(3)=tgf val(4)=tg call printv(imsg,val,zval,4) if (o_corr) then write(imsg,1010) else write(imsg,1000)ipstep endif call vec_dcopy(point,sw,nn) c BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB c B dstep is a Branin-Step, if o_nn is true !!! B c BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB if (o_nn) then bn=bnorm(adj,gradient,gmat,dstep,nn) if(o_even) then call vec_dmult_constant(dstep,nn,-pstepl,dstep) else call vec_dmult_constant(dstep,nn,pstepl,dstep) endif endif c ccccccccccccccccccc call vec_dmult_add(point,dstep,nn,1.0d0,point) call vec_dcopy(dstep,cold,nn) eold=energy irgf_para(1)=istep irgf_para(2)=ipstep irgf_para(3)=itmin irgf_para(4)=itmax call recreate(isto,z_store) open(isto,file=z_store) write(isto,*)fstep write(isto,*)proj write(isto,*)told write(isto,*)t1old write(isto,*)cold write(isto,*)xold write(isto,*)gold write(isto,*)hold write(isto,*)eold write(isto,nml=store) close(ipnt,status='delete') open(ipnt,file=z_FoP,status='new') write(ipnt,*)point if (istep.eq.1) call print_switch(ires) write(ires,1040)istep,itmax write(imsg,1040)istep,itmax if (o_bif) write(ires,1020) if (o_turn) write(ires,1030) if (o_corr) then write(ires,1010) else write(ires,1000)ipstep endif write(ires,1060)dmnorm(dstep,gmat,nn) if (dmnorm(dstep,gmat,nn).gt.pstepl*1.05d0) then write(imsg,*)'Oversized' endif do i=1,nn write(ires,1070)sw(i),dstep(i) enddo write(ires,1100)(isig(i),i=1,3) write(ires,1050)gnorm,dnst,crit_stop,rgnorm,crit_pc write(imsg,1050)gnorm,dnst,crit_stop,rgnorm,crit_pc call closext call closys if (istep.eq.itmax) then irgfmain=10 return endif if ((dnst.le.crit_stop).and.(ipstep.ge.itmin)) then irgfmain=10 return endif irgfmain=2 if (o_corr) then irgfmain=irgfmain+4 else irgfmain=irgfmain+2 endif if (o_req_ex) irgfmain=irgfmain+1 return 1000 format(' predictorstep',i4) 1010 format(' correctorstep') 1020 format(' bifurcation criterion satisfied') 1030 format(' turning point passed') 1040 format(1X,20('='),/,' step:',i4,'/',i4,/,1x,20('=')) 1050 format(/36x,'Value',9x,'Criterion', c /1x,'Length of the gradient: ',f14.6,13x,'[0.0]', c /1x,'Stopping criterion: ',f14.6,4x,f14.6, c /1x,'Predictor-Corrector Crit.:',f14.6,4x,f14.6) 1060 format(1x,'Length of Step',f9.6, c /9x,'Point Step') 1070 format(4x,f12.7,' ',f10.7) 1080 format(1x,'Euclidian metric') 1090 format(1x,'Non-Euclidian metric') 1100 format(1x,'signature of the Hessian: ',i2,'+ ',i2,'- ',i2,'=') 1110 format(1x,'Binary') 1120 format(1x,'Messages') 1130 format(1x,'Exact Hessian') 1140 format(1x,'S Norm ',f9.6,'\t',A20) 1150 format(1x,'Rayleigh ',f9.6,'\t',A20) end