      PROGRAM A22Vri3         
! Programm-Teil 3: organize and execute the VRI search
! Programm bestimmt Suchrichtung der Newtontrajektorie, die den
! VRI-Punkt enthält
!
! W.Quapp   2009
! 
!      n = Dimension
! length = Kettenlaenge
!    eps = Genauigkeit des RGF-Verfahrens
!  kette = Kette zwischen den Punkten xstart und xend=vri
! lkette = chain from part 2, length le
!    dir = Suchrichtung des RGF-Verfahrens = rsel
!      
      INTEGER nn,le,N3,I5
      DOUBLE PRECISION eps
      PARAMETER (nn=58,N3=22,le=12,I5=3)
!
      DOUBLE PRECISION lkette(le+1,nn),kette(le+1,nn),
     1  dir(nn),xstart(nn),xend(nn),vri(nn),
     2  vridir(nn),Agminx,AgMin,xx,
     3  grad(nn),H(nn,nn),Adj(nn,nn),Ag(nn),
     4  vriold(nn),dmnorm,norm,tmin,
     5  gradj(I5+1,nn),Hj(I5+1,nn,nn),
     6  ginvj(I5+1,nn,nn),gmatj(I5+1,nn,nn),
     7  ginv(nn,nn),gmat(nn,nn),tgt,sgn
      INTEGER indexJ,istatus,iter,j1,k1,le22,jj,icalc
ccccccccccccccccccccccccccccccccccccccccc
      OPEN(7,FILE='chain.dat',status='unknown',access='append')
      OPEN(8,FILE='a22Vkette.dat')
      OPEN(9,FILE='a22VketNrj.dat')
      OPEN(10,FILE='a22Vket3.dat')
      OPEN(11,FILE='a22Vdir.dat')
      OPEN(14,FILE='a22Vgrad.dat')
      OPEN(15,FILE='a22Vhesse.dat')
      OPEN(16,FILE='a22Veps.dat')
      OPEN(17,FILE='a22Viter.dat')
      OPEN(18,FILE='a22Vistat.dat')
      OPEN(19,FILE='a22Vicalc.dat')
      OPEN(28,FILE='a22VketNrj1.dat')
      OPEN(30,FILE='a22Vvri.dat')
      OPEN(33,FILE='a22Vmin.dat')
      OPEN(48,FILE='a22Vgmatr.dat')
c protocol file: all interesting output of the vri-search
      OPEN(44,FILE='a22Vprotocol.txt',status='unknown',access='append')
cccccccccccccccccccccccccccccccccccccccc
c
       rewind 16
       read(16,*) eps
       j=1
       rewind 9
       write(9,*) j  
       rewind 28
       read(28,*) j1,k1
       rewind 19
       read(19,*) icalc,sgn
       rewind 18       
       read(18,*) istatus
       rewind 17
       read(17,*) iter
       rewind 33
       read(33,*) AgMin, Agminx, xx
c start vri:
       if(icalc.eq.0) then
        icalc=2
        rewind 19
        write(19,*) icalc,sgn
        goto 77
       endif
c loop  vri:
       if(icalc.eq.1) then
        icalc=2
        rewind 19
        write(19,*) icalc,sgn
        goto 77
       else
        icalc=1
        rewind 19
        write(19,*) icalc,sgn
        goto 88
       endif 
ccccc       
77     continue
ccccc
       j1=j1+1
       if(j1.eq.le+2) then
         istatus=2
         rewind 18
         write(18,*) istatus
         goto 55
       else 
         istatus=4
         rewind 18
         write(18,*) istatus
       endif
       
       write(*,*) '              j1 in  iteration=  (',j1,iter,')'
       write(44,*) '             j1 in  iteration=  (',j1,iter,')'
       rewind 28
       write(28,*) j1,k1
       rewind 8
       Do 23 k=1,le+1
        read(8,*) (lkette(k,i),i=1,nn)
c         Do 22 i=1,nn 
c22       rkette(k,i)=lkette(k,i)
23     continue 
         rewind 30
         read(30,*) (vri(i),i=1,nn)
         read(30,*) (vriold(i),i=1,nn)
50    FORMAT(1X,3F20.10)         
        if(j1.eq.2) then
          WRITE(*,*) ' vri '
          WRITE(*,50) (vri(i),i=1,5)
          WRITE(*,50) (vri(i),i=6,nn)
        endif
c use only xstart on chain, xend will be the former vri point
         DO 10 i=1,nn
           xstart(i)= lkette(j1,i)
cc         xend(i)  = lkette(le+1,i)
10       CONTINUE
         DO 173 i=1,nn
         xend(i)  =vriold(i)
173      CONTINUE
2      FORMAT(3F20.14)
       CALL straight_chain(kette,xstart,xend,nn,le)
c staerke Kettendichte am Ende, um bei VRI genauer zu werden
         DO 174 i=1,nn
          xstart(i)  = kette(le-2,i)
174      CONTINUE
       CALL straight_chain(kette,xstart,xend,nn,le)           
c Use only last I5+1 nodes of the chain up to VRI
C      I5=3  is put above 
c Muss in Teil 4 analog eingestellt sein !!
      rewind 10
      DO 313 k=le-I5,le+1
        WRITE(10,2) (kette(k,i),i=1,nn)
313   CONTINUE
      goto 55  
cccccccccccccccc
88    continue
cccccccccccccccc
       rewind 30
       read(30,*) (vri(i),i=1,nn)
       read(30,*) (vriold(i),i=1,nn)
c  I5 siehe oben !!
       rewind 10
       DO 32 jj=1,I5+1
         READ(10,*) (kette(jj,i),i=1,nn)
32     CONTINUE
       rewind 14
       DO 111 jj=1,I5+1
111      read (14,*) (gradj(jj,i),i=1,nn)
       rewind 15
       DO 110 jj=1,I5+1
         Do 113 k=1,nn
113        read (15,*) (Hj(jj,k,i),i=1,nn)
110    continue
       rewind 48
       DO 114 jj=1,I5+1
         Do 766, k=1,nn
766         READ(48,*) (ginvj(jj,k,I),I=1,nn)
            READ(48,*)
            Do 768, k=1,nn
768         READ(48,*) (gmatj(jj,k,I),I=1,nn)
            READ(48,*)
114    continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      indexJ=0  
      DO 101 jj=1,I5 + 1 
c  use current values of grad, metric, Hessian 
       DO 115 i=1,nn
115       grad(i)= gradj(jj,i)
       Do 116 k=1,nn
       DO 116 i=1,nn 
          gmat(k,i)= gmatj(jj,k,i)
          ginv(k,i)= ginvj(jj,k,i)
          H(k,i)   =    Hj(jj,k,i)
116    continue
c shift Hessian !!
       call vec_dmult_constant(H,nn*nn,4.0d0,H)       
       tgt=dmnorm(grad,ginv,nn)
       write(*,*) '|grad|=   ',tgt
       write(44,*)'|grad|=   ',tgt
c test to avoid points to near to SP
c        IF(tgt.gt.0.02d0) THEN      ! avoid SP use for VRI 
         CALL mat_adj(H,nn,Adj)
         CALL matmult(Adj,nn,grad,1,Ag,nn)
c H, g are in bohr, radian ( kette is in angstroem, degree )
C Metric norm: , Ag is 1 x contravariant:
         tmin=dmnorm(Ag,gmat,nn) 
         write(44,*) '  tmin with      gmat-norm               ',tmin
         write(*,*)  '  tmin with      gmat-norm               ',tmin
         IF(tmin .LT. AgMin ) THEN
           AgMin=tmin
           rewind 33
           write(33,2) tmin, AgMinx, xx
           write(*,*)'          new Ag  -  minimum  = ',tmin 
           write(44,*)'         new Ag  -  minimum  = ',tmin
           do 102 i=1,nn
               vri(i)=kette(jj,i)
102            vridir(i)=gradj(jj,i)
           indexJ=jj  
           ENDIF
c        ENDIF
101    CONTINUE
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       if( (indexJ.ne.0) ) then
         write(*,*) 'vri   '
         write(*,50) (vri(i),i=1,5)
         write(*,50) (vri(i),i=6,nn)
         WRITE(*,*) 'vriold'
         write(*,50) (vriold(i),i=1,5)
         WRITE(*,50) (vriold(i),i=6,nn)
         write(*,*) 'vridir' 
         write(*,50) (vridir(i),i=1,5)
         write(*,50) (vridir(i),i=6,nn)
ccc          
         write(44,*) 'vri   '
         write(44,50) (vri(i),i=1,5)
         write(44,50) (vri(i),i=6,nn)
         WRITE(44,*) 'vriold'
         write(44,50) (vriold(i),i=1,5)
         WRITE(44,50) (vriold(i),i=6,nn)
         write(44,*) 'vridir'
         write(44,50) (vridir(i),i=1,5)
         write(44,50) (vridir(i),i=6,nn)
ccc
         rewind 30
         write(30,2) (vri(i),i=1,nn)
         write(30,2) (vriold(i),i=1,nn)
         rewind 11
         write(11,2) (vridir(i),i=1,nn)
ccc
         WRITE(7,*) '(* vri-point (phi,psi), iteration (',j1,iter,')'
         WRITE(7,*) '{',vri(20),',',vri(44),' } *)'
       endif
c end jj cycle for any j1  !
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
55    CONTINUE
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       rewind 18
       write(18,*) istatus
       call exit(istatus)
      END
!
!  Kette der Länge le zwischen Punkten x und y
!        
      SUBROUTINE straight_chain(kette,x,y,nn,le)
      DOUBLE PRECISION kette(le+1,nn),x(nn),y(nn)
      DO 20 j=1,le+1
        DO 10 i=1,nn
                kette(j,i)=x(i)+((j-1)*(y(i)-x(i)))/le 
10      CONTINUE
20    CONTINUE
      END
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*8 function dmnorm(vec,gm,nn)
c   metric norm of VEC with resp. to metric GM
      real*8 vec,vec_mdot,gm
      integer nn
      dimension vec(nn),gm(nn,nn)
          dmnorm=dsqrt(vec_mdot(vec,gm,vec,nn))
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
! Berechnet adjunkte Matrix A zu M, in 3D only ! 
      subroutine mat_adj(matin,ndim,matout)
c adjoint matrix of matin
      real*8 matin,matout,matscr,det
      integer ndim,i,j,k,l,ii,jj
      dimension matin(ndim,ndim),matout(ndim,ndim),
     +  matscr(ndim-1,ndim-1)
      call vec_dinit(matscr,(ndim-1)*(ndim-1),0.0d0)
      do 10, i=1 , ndim
        do 20, j=1 , ndim
             do 30, ii=1 , ndim-1
                do 40, jj=1 , ndim-1
                        if (ii.lt.i) then
                         k=ii
                        else
                         k=ii+1
                        endif
                        if (jj.lt.j) then
                         l=jj
                        else
                         l=jj+1
                        endif
                        matscr(ii,jj)=matin(k,l)
                       if(dabs(matscr(ii,jj)).lt.1.0d-63) then
                         matscr(ii,jj)= 1.0d-63
                       endif
 40               continue
 30             continue
                matout(i,j)=det(matscr,ndim-1)
                if(dabs(matout(i,j)).lt.1.0d-63) goto 20
                matout(i,j)= (-1.0d0)**(i+j)*matout(i,j)
 20     continue
 10   continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccc
       real*8 function det(matein,nn)
c   calculation of the determinant of MATEIN
       integer i,nn
       real*8 matein,mataus,w
       dimension matein(nn,nn),mataus(nn,nn),w(nn)
       call householder(matein,mataus,w,nn,nn)
       det=1.0d0
       do 10,i=1,nn
       if (dabs(mataus(i,i)).lt.1.0d-63) then
          det= 1.0d-63
          goto 20
       endif
       if (dabs(det) .lt.1.0d-63) then
          det= 1.0d-63
          goto 20
       endif
 10    det=det*mataus(i,i)
       if (iand(nn,1).eq.0) det=-det
 20    continue
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine householder(matein,mataus,w,n,sp)
c    householder transformation
       integer k,l,m,n,sp
       real*8 alpha,rho,sk,s,w,matein,mataus
       dimension matein(n,sp),mataus(n,sp),w(n)
       call vec_dcopy(matein,mataus,n*sp)
       do 10,k=1,n-1
        s=0.0d0
        do 20,l=k,n
         s=s+mataus(l,k)*mataus(l,k)
 20     continue
        if (mataus(k,k).lt.0) then
         alpha=dsqrt(s)
        else
         alpha=-dsqrt(s)
        endif
        rho=dsqrt(s-mataus(k,k)*mataus(k,k)+
     *       (mataus(k,k)-alpha)*(mataus(k,k)-alpha))
        if(rho .lt. 1.0d-63) then
         rho= 1.0d-63
            w(k)= 0.0d0
            goto 22
        endif
        w(k)=(mataus(k,k)-alpha)/rho
 22     continue
        do 30,l=k+1,n
         w(l)=mataus(l,k)/rho
 30     continue
        do 40,l=k,sp
         sk=0.0d0
         do 55,m=k,n
          sk=sk+w(m)*mataus(m,l)
 55      continue
         do 60,m=k,n
           mataus(m,l)=mataus(m,l)-2*sk*w(m)
 60      continue
 40     continue
 10    continue
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccc
!
! Gibt die Norm von v aus
!
      DOUBLE PRECISION FUNCTION norm(v,nn)
      DOUBLE PRECISION v(nn),s
      s=0.d0
      DO 10 j=1,nn
        s=s+v(j)**2
10    CONTINUE
      IF (s.LT.1.D-30) THEN
        norm=0.d0
      ELSE
        norm=DSQRT(s)
      ENDIF
      END
!
! Gibt die Distanz zweier Punkte v und w aus
!
      DOUBLE PRECISION FUNCTION distance(v,w,nn)
      DOUBLE PRECISION v(nn),w(nn),s      
      s=0.d0
      DO 10 j=1,nn
        s=s+(v(j)-w(j))**2
10    CONTINUE
      distance=DSQRT(s)
      END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine tri_mat(dm1,nr1,nc1,dm2,nc2,dm3,nc3,dres)
c     dm1[nr1,nc1]*dm2[nc1,nc2]*dm3[nc2,nc3]=dres[nr1,nc3]
      real*8 dm1,dm2,dm3,dres,r1
      integer nr1,nc1,nc2,nc3
      dimension dm1(nr1,nc1),dm2(nc1,nc2),dm3(nc2,nc3),
     +          r1(nr1,nc2),dres(nr1,nc3)
      call matmult(dm1,nr1,dm2,nc2,r1,nc1)
      call matmult(r1,nr1,dm3,nc3,dres,nc2)
      return
      end
ccccccccccccccc
       subroutine mat_dcopy(arra1,arra2,ndim,mdim)
c  Copies matrix ARRA1 to matrix ARRA2
       real*8  arra1,arra2
       integer i,j,ndim,mdim
       dimension arra1(ndim,mdim),arra2(ndim,mdim)
       do 10,i = 1,ndim
       do 20,j = 1,mdim
        arra2(i,j) = arra1(i,j)
 20    continue
 10    continue
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
      real*8 function vec_mmdot
     &(nri,matin,din,nra,matout,dout,dim,dmet)
c     multiply the (nri) col of (matin) and the (nra) col of
c     (matout). (dim) is the number of the rows with metric dmet
      real*8 matin,matout,prod,v1,v2,dmet,vec_mdot
      integer nri,din,nra,dout,dim
      dimension matin(dim,din),matout(dim,dout),dmet(dim,dim),
     +   v1(dim),v2(dim)
ccc      if ((nri.gt.din).or.(nra.gt.dout)) call cerr(-100)
      call vec_mat_centry(v1,matin,nri,din,dim,'o')
      call vec_mat_centry(v2,matout,nra,dout,dim,'o')
      prod=vec_mdot(v1,dmet,v2,dim)
      vec_mmdot=prod
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine vec_dcopy(arra1,arra2,ndim)
c  Copies ARRA1 to ARRA2
      real*8 arra1,arra2
      integer i,ndim
      dimension arra1(ndim),arra2(ndim)
      if(ndim .le. 0) return
      do 10 i = 1,ndim
           arra2(i) = arra1(i)
 10   continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine vec_mat_centry(vec,mat,npos,ndim,ncol,zkind)
c Copies VEC 'i'n or 'o'ut the NPOS column of matrix MAT
        real*8 vec,mat
        integer npos,ndim,ncol,i
        character zkind
        dimension vec(ndim),mat(ndim,ncol)
        if(npos.gt.ncol) return
        if (zkind.eq.'i') then
         do i=1,ndim
           mat(i,npos)=vec(i)
         end do
        endif
        if (zkind.eq.'o') then
         do i=1,ndim
           vec(i)=mat(i,npos)
         enddo
        endif
        return
        end
cccccccccccccccccccccccccccccccccccccccccccccccccccc
       real*8 function vec_mdot(vec1,metric,vec2,dim)
c  metric scalar product of VEC1 and VEC2 with respect to METRIC
       integer dim
       real*8 vec1,metric,vec2,cov2,vec_ddot
       dimension vec1(dim),vec2(dim),cov2(dim),metric(dim,dim)
       call matmult(metric,dim,vec2,1,cov2,dim)
       vec_mdot=vec_ddot(vec1,cov2,dim)
       return
       end
ccccccccccccccccccccccccccccccccccccc
       real*8 function vec_ddot(arra1,arra2,ndim)
c   cartesian scalar product of ARRA1 and ARRA2
       real*8  arra1,arra2
       integer ndim,i
       dimension arra1(ndim),arra2(ndim)
       vec_ddot = 0
       if(ndim .le. 0)return
       do 70 i = 1,ndim
         vec_ddot = vec_ddot + arra1(i) * arra2(i)
  70   continue
       return
       end
ccccccccccccccccccccccccccccccccccccccc
       subroutine matmult(m1,d1,m2,d2,mout,dim)
c  matrix multiplication:  m1[d1,dim]*m2[dim,d2]=mout[d1,d2]
       real*8 m1,m2,mout,swap
       integer d1,d2,dim,i,j,k
       dimension m1(d1,dim),m2(dim,d2),mout(d1,d2),swap(d1,d2)
       call vec_dinit(swap,d1*d2,0.0d0)
       do 10,i=1,d1
        do 20,j=1,d2
         do 30,k=1,dim
           swap(i,j)=swap(i,j)+m1(i,k)*m2(k,j)
 30      continue
 20     continue
 10    continue
       call vec_dcopy(swap,mout,d1*d2)
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine vec_dinit(arra1,ndim,wert)
c Set all elements of ARRA1 to the double value of WERT
       real*8   arra1,wert
       integer  ndim,i
       dimension arra1(ndim)
       if(ndim .le. 0) return
       do i=1,ndim
         arra1(i)=wert
       end do
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
      subroutine mat_trans(matin,nrow,ncol,matout)
c     transpose MATOUT of matrix MATIN
      real*8 matin,matout,swap
      integer nrow,ncol,i,j
      dimension matin(nrow,ncol),matout(ncol,nrow),swap(ncol,nrow)
      do 10,i=1,nrow
        do 20,j=1,ncol
                swap(j,i)=matin(i,j)
   20   continue
   10 continue
      call mat_dcopy(swap,matout,ncol,nrow)
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine vec_mrowcopy(nri,matin,din,nra,matout,dout,dim)
c  copies the (nri) row of (matin) with (din) row and (dim) col
c  to the (nra) row of (matout) with (dout) row and (dim) col
      real*8 matin,matout
        integer nri,din,nra,dout,dim,i
        dimension matin(din,dim),matout(dout,dim)
ccc      if ((nri.gt.din).or.(nra.gt.dout)) call cerr(17)
      do 10,i=1,dim
        matout(nra,i)=matin(nri,i)
   10 continue
        return
        end
cccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine vec_madd
     * (nr1,mat1,d1,fact1,nr2,mat2,d2,fact2,nrx,matx,dx,dim)
c     matx(nrx)=fact1*mat1(nr1)+fact2*mat2(nr2)
c     d1,d2,dx - columns
c     dim      - rows
       real*8 mat1,mat2,matx,fact1,fact2,vec
       integer nr1,nr2,nrx,d1,d2,dx,dim,i
       dimension mat1(dim,d1),mat2(dim,d2),matx(dim,dx),vec(dim,1)
       call vec_dinit(vec,dim,0.0d0)
       do 10,i=1,dim
          vec(i,1)=fact1*mat1(i,nr1)+fact2*mat2(i,nr2)
 10    continue
       call vec_mat_centry(vec,matx,nrx,dim,dx,'i')
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine vec_dmult_constant(arra1,ndim,zahl,arra2)
c   ARRA2 = ZAHL * ARRA1
      real*8  arra1,arra2,zahl
      integer ndim,i
      dimension arra1(ndim),arra2(ndim)
      if(ndim .le. 0)return
      do i = 1,ndim
         arra2(i) = arra1(i)*zahl
      end do
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine mat_diag(d,n,v)
c create a NxN diagonal matrix D with entries V
       real*8  d,v
       integer n,i
       dimension d(n,n)
              if(n .le. 0) return
               call vec_dinit(d,n*n,0.0d0)
               do i=1,n
                 d(i,i)=v
               end do
               return
              end
ccccccccccccccccc