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  
