c c TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT c T mTASC - Version 1 August 31, 2001 by Wolfgang Quapp T c T (see TCA 105 (2000) p.145-155 ) T c TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT c c based on RGF modular version c Revision 2.2, Aug 31, 2001 c by c Michael Hirsch / Wolfgang Quapp c Uni. Leipzig c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 TASC real*8 crit_stop,crit_pc,pstepl,some_real,eold integer nn,irgfmain,istat,itmin,itmax,irgf_para,narg,iargc integer job_id,job integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene integer 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 job integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene real*8 crit_stop,crit_pc,dsqrt,hessian,dold,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 Cstep,SLS,energy,eold c real*8 some_real,ray,snorm,dray,dsn 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 rsel(nn),hessian(nn,nn),gradient(nn),fstep(nn) dimension 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),dold(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) c 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,o_SL,igstep,eold,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,*)dold 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 igstep=0 o_SL=.false. 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) read(iwil,*,iostat=iochk)bmat read(iene,*,iostat=iochk)energy c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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) if(.not.o_SL) then call err(i_dfp(hold,dg,dx,nn,hessian),'DFP',3) call message('(1x,''=== updated Hessian ==='')',30) endif endif c 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) if (o_msg) write(ires,*) ' eigenvalues: ' if (o_msg) write(ires,*) (eival(i),i=1,nn) 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 dnst=newton(hessian,gradient,dstep,gmat,nn) c 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) c if (o_even) 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 c c BigLLLLLLLOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOPPPPPPPPPP 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 c GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG if(o_even .and. (.not.o_SL) .and. eival(1).lt.0 c .and. Dabs(eival(1)).gt. eival(2) ) then c Do Steepest Descent o_corr=.true. SLS=1.0d0 if(gnorm .lt. pstepl) SLS=pstepl / gnorm if(gnorm .gt. 2.5d0*pstepl) SLS=2.5d0*pstepl / gnorm o_corr=.true. igstep=igstep+1 istep=istep-1 write(imsg,*)' No ',igstep,' of Steepest Descent ' call vec_dmult_constant(gradient,nn,-SLS,dstep) cc End Gradient search else c GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG c 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 rgnorm=dmnorm(sw,gmat,nn) if (.not.o_SL) then call gettang(rgeq,gmat,tang,nn) else cc use old tangent for SL search call vec_dcopy(told,tang,nn) endif c call vec_dinit(sw,nn,0.0d0) 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) c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c if ((rgnorm .gt. crit_pc) .and. (.not.o_SL) ) then c c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c Corrector call linsolve(dk,dstep,nn) Cstep=dmnorm(dstep,gmat,nn) if(Cstep.gt.pstepl) then call vec_dmult_constant(dstep,nn,pstepl/Cstep,dstep) write(ires,*) '! Reduced Corrector ' endif if (o_recorr) c call recorrect(dstep,dold,gmat,o_corr,nn) if (o_corr) call message('(1x,''recorrect'')',16) o_corr=.true. o_SL=.false. c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC else c PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPP c Predictor SLS=pstepl if(.not.o_SL) then o_SL=.true. ipstep=ipstep+1 else o_SL=.false. istep=istep-1 call search_line(o_even,gold,gradient,eold,energy, c xold,point,eival(1),SLS,tang,nn) if(SLS.gt.pstepl) SLS= pstepl SLS= SLS / dmnorm(tang,gmat,nn) endif call vec_dmult_constant(tang,nn,SLS,dstep) cc End TASC SL_search c LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL c 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 c TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT C make TASC - SEARCH - DIRECTION if( (.not.o_SL) .and. (.not.o_corr) .and. (ipstep.gt.0) c .and. Dabs(eival(1)).lt.0.9d0*eival(2) ) then call ortcompl(tang,proj,unit,nn) write(ires,*) ' TASC - step ' endif if ( Dabs(eival(1)) .gt. 0.9d0*eival(2) ) then write(ires,*) 'RGF - step to Don Quixote SP' o_SL=.false. endif c TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT endif c of SteepestDescent - Predictor / Corrector - Loop endif endif val(1)=detk val(2)=tgt val(3)=tgf val(4)=tg call printv(imsg,val,zval,4) if (o_corr) then if ( istep .eq.1) then write(imsg,*) ' Steepest Descent - step ', igstep else write(imsg,1010) endif else write(imsg,1000)ipstep endif c c Last Step: if( (dnst.le.crit_stop) .and. (ipstep.ge.itmin) ) then write(ires,*) '==================== ' write(ires,*) 'Do last step by Newton-Raphsen ' dnst=newton(hessian,gradient,dstep,gmat,nn) call vec_dcopy(point,xold,nn) call vec_dmult_add(point,dstep,nn,1.0d0,point) o_SL=.false. c Last Step else if (o_corr.or.(istep.eq.1)) then call vec_dcopy(point,xold,nn) call vec_dmult_add(point,dstep,nn,1.0d0,point) else if(o_SL) then call vec_dcopy(point,xold,nn) call vec_dmult_add(point,dstep,nn,1.0d0,point) else call vec_dmult_add(xold,dstep,nn,1.0d0,point) endif endif endif c eold=energy call vec_dcopy(dstep,dold,nn) call vec_dcopy(gradient,gold,nn) c 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,*)dold 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 (ipstep.eq.0) 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( (dnst.le.crit_stop) .and. (ipstep.ge.itmin) ) then write(ires,*) ' Newton-Step ' else if (o_corr) then if ( istep .eq.1) then write(ires,*) ' Steepest Descent - step ', igstep else write(ires,1010) endif else write(ires,1000)ipstep endif endif write(ires,1060)dmnorm(dstep,gmat,nn) if (dmnorm(dstep,gmat,nn).gt.pstepl*1.05d0) then write(imsg,*)'Oversized' endif if (.not.o_SL .or. o_corr) then do i=1,nn write(ires,1070)xold(i),dstep(i) enddo endif 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 if ((dnst.le.crit_stop).and.(ipstep.ge.itmin) ) then write(ires,*) ' Last Point by Newton-Raphsen' do i=1,nn write(ires,1070) point(i) enddo write(ires,*) ' Program has stopped ' endif call closext call closys if (istep.eq.itmax) then irgfmain=10 return endif if ((dnst.le.crit_stop).and.(ipstep.ge.itmin) c .and. ((.not.o_SL ).or. o_corr)) 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(' RGF-bifurcation criterion satisfied') 1030 format(' RGF-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