        PROGRAM A22Vri2
!
!  Part 2 :   rgf to dir=rsel  
! 
!   W.Quapp, Dezember 2009
!
!    CCC for Alanine potential CCCCCCCCCCCCCCCCCC
! energy + gradient + Hessian + metric from GamessUS
! Programm bestimmt Suchrichtung der Newtontrajektorie, die den
! VRI-Punkt enthält
!
!     nn = Dimension,  N3 = number of atoms
!     le = Chainlength (even integer)
!   pstl = pred-Steplength
!    eps = Convergence in RGF corrector,
!                      is useless for predictors only
!  chain = Chain between points xstart and xend,
!                      xend is useless for predictors only
!    dir = Searchdirection RGF method = rsel in RGF corrector
!
! scriHCNntVRI: is the working code in Linux, combining the
! followings FORTRAN programms, which use following steps
!
! (1) Start: Initial settings, with
!  (i) Fix start- and end point of the chain
!  (ii) Fix searchdirection "dir", and initial guess of "VRI"
!  (iii) Put straigt chain of "length" between both points
!    This parts:
! (2) To "dir" calculate corresponding NT
!  use K-mimic; other version: use Ag for a predictor step

! (3) Calculate VRI approximation,
!     use its gradient for new "dir" in (2)
! (4) Organize points for Gamess-US input for (2) or (3)
! A22vriReadOUT: translates output into usable input files
! (5) Organize Iterations up to convergence of VRI
!     by repeating of (2)-(4)
! 
      INTEGER nn,N3,le,iterNTmax
      DOUBLE PRECISION pst
      PARAMETER (nn=58,N3=22,le=12,iterNTmax=6,pst=0.2d0)
!
      DOUBLE PRECISION kette(le+1,nn),rkette(le+1,nn),eps,tgt,tt,pla,
     1  rsel(nn),rc(nn),grad(nn),H(nn,nn),Proj(nn,nn),Ag(nn),tmin,
     2  xstart(nn),xend(nn),sw(nn),prestep(nn),Adj(nn,nn),Agmin,
     3  norm,gmat(nn,nn),dstep(nn),detk,tminNT,AgmiNT,factor,
     4  projNm1(nn-1,nn),dk(nn,nn+1),rgeq(nn-1,nn),Help(nn,nn),
     5  rgrad(nn-1),tcov(nn),tang(nn),told(nn),
     6  unit(nn,nn),dmnorm,vec_mdot,vec_ddot,det,pstl,
     7  t1old(nn),t1(nn),rgnorm,point(nn),eival(nn),eivec(nn,nn),
     8  gradcont(nn),ginv(nn,nn),energy,sgn,detHe
      INTEGER k,j,i,jnr,istatus,icalc,iter,iterNT,Natoms(N3)
      logical omat_inv
      Character*2  CharAtoms(N3)
      Character*5  Zcoor(nn+2)
      Character*70 zmatLines(35)
      Data Zcoor/'ch2  ','cc3  ','cch3 ','oc4  ','occ4 ','dih4 ',
     1   'hc5  ','hcc5 ','dih5 ','hc6  ','hcc6 ','dih6 ','nc7  ',
     2   'ncc7 ','dih7 ','cn8  ','cnc8 ','dih8 ','cc9  ','ccn9 ',
     3   'dih9 ','oc10 ','occ10','dih10','hc11 ','hcn11','dih11',
     4   'hn12 ','hnc12','dih12','cc13 ','ccn13','dih13','hc14 ',
     5   'hcc14','dih14','hc15 ','hcc15','dih15','hc16 ','hcc16',
     6   'dih16','nc17 ','ncc17','dih17','cn18 ','cnc18','dih18',
     7   'hc19 ','hcn19','dih19','hn20 ','hnc20','dih20','hc21 ',
     8   'hcn21','dih21','hc22 ','hcn22','dih22'/
      Data CharAtoms/'h','c','c','o','h','h','n','c','c','o','h','h',
     1               'c','h','h','h','n','c','h','h','h','h'/
      Data Natoms/1,6,6,8,1,1,7,6,6,8,1,1,6,1,1,1,7,6,1,1,1,1/
c     pure and applied chemistry, 51, 1 (1979):
c     angstroms per bohr
       data toang  /0.52917706d0/
ccccccccccccccccccccccccccccccccccccccccc      
      OPEN(8,FILE='a22Vkette.dat')
      OPEN(7,FILE='chain.dat',status='unknown',access='append')
      OPEN(9,FILE='a22VketNrj.dat')
      OPEN(11,FILE='a22Vdir.dat')
      OPEN(13,FILE='a22Vener.dat')
      OPEN(53,FILE='a22VEall.dat',status='unknown',access='append')
      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(25,FILE='a22Vtold.dat')
      OPEN(26,FILE='a22Vt1old.dat')
      OPEN(28,FILE='a22VketNrj1.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
2     FORMAT(3F20.14)
      N=nn
      j1=1
      k1=1  
c     j1 for iteration in part Q3, k1 is obsolet
      rewind 28
      write(28,*) j1,k1 
      rewind 16
      read(16,*) eps
c  eps is criterion for rgf convergence to given rsel 
       rewind 9
       Read(9,*) jnr
       rewind 19
       Read(19,*) icalc,sgn
c jnr is the currend node number
      If( (jnr.eq.1) .and. (icalc.eq.0) ) then
         write(*,*)
         write(*,*) 'eps = ',eps
         write(44,*)'eps = ',eps
       endif
       rewind 8
       DO 11 k=1,le+1
          Read(8,*) (kette(k,i),i=1,nn)
11      CONTINUE
       rewind 17    
       Read(17,*) iter,iterNT
       if(iter.eq.0) then
          pstl=pst
       else
          pstl=1.755555d0*pst
       endif 
      if(iterNT.ge.iterNTmax) then
        write(*,*)
        write(*,*)'Stop: NT-node not convergent at jnr',jnr
     1           ,'  in Iteration',iter 
        write(*,*)   
        write(44,*)
        write(44,*)'Stop: NT-node not convergent at jnr',jnr
     1           ,'  in Iteration',iter 
        write(44,*)
        if(jnr.lt.le+1) then 
           jnr=jnr+1
         else 
           istatus=0
           rewind 18
           write(18,*) istatus
           sgn=-1.0d0
           icalc=0
           rewind 19
           write(19,*) icalc,sgn
           goto 66           
        endif
        iterNT=0
        write(17,*) iter,iterNT
       endif

       if(icalc.eq.1) then 
        write(44,*) 'kette(',jnr,') angstroem,degree in Vri2 '
        Do 125, i=1,5 
125      write(44,*) i,'    ',Zcoor(i), kette(jnr,i)
        Do 12, i=6,nn
12       write(44,*) i,'   ',Zcoor(i+1), kette(jnr,i)
       endif 
c c ccccccccccccccccccccccccccccccccccccccccccccccccc
c (bohr,radian)<->(angstroem,degree) for internal coordinates
c  current point  (angstroem,degree) ->(bohr,radian)
        torad=dacos(-1.0d0)/180.0d0
        btoa =0.52917706d0
        point(1)=kette(jnr,1)/btoa
        point(2)=kette(jnr,2)/btoa
        point(3)=kette(jnr,3)*torad
        point(4)=kette(jnr,4)/btoa
        point(5)=kette(jnr,5)*torad
       DO 13 k=2,18
         point(3*k  )= kette(jnr,3*k  )/btoa
         point(3*k+1)= kette(jnr,3*k+1)*torad
         point(3*k+2)= kette(jnr,3*k+2)*torad
13      continue
         point(57)=kette(jnr,57)/btoa
         point(58)=kette(jnr,58)*torad
c
c in corrector loop with use of gradient and hessian, there
c are used bohr and radian. For input of Gamess, however, it is
c angstroem and degree !!!  
c     write(44,*)'point in bohr,radian, in *Vri2.f '
c     write(44,2)(point(i),I=1,nn)
c c ccccccccccccccccccccccccccccccccccccccccccccccccc
      call mat_diag(unit,nn,1.0d0)
c
      rewind 11
      Read(11,*) (rsel(i),i=1,nn)
      if(jnr.eq.1) then
c         rewind 14
c         if(iter.lt.1) then 
c           Read(14,*) (rsel(i),i=1,nn)  
c           write(44,*) 'r_selected in vri2 test: bei jnr 1 = gradient '
c         else
           write(44,*) 'r_selected = grad_vri  '
c         endif  
         write(44,2)(rsel(i),i=1,5)
         write(44,2)(rsel(i),i=6,nn)
      endif
c      
      rewind 18 
      read(18,*) istatus
      if(istatus.eq.10) then
         write(44,*)'Abbruch in vri2 cycle'
         write(*,*) 'Abbruch in vri2 cycle'
         goto 66 
      else
         istatus=2
         rewind 18
         write(18,*) istatus
      endif 
c
       rewind 25
       Read(25,*) (told(i),i=1,nn)
       rewind 26
       Read(26,*) (t1old(i),i=1,nn)
c              
       if(icalc.eq.0) then
        icalc=1
        write(19,*) icalc
         goto 66
       else
        icalc=0
        write(19,*) icalc
        rewind 13
        Read(13,*) energy
        write(44,*)'Energy ',energy,' at iter ',iter,
     1    'kete(',jnr,')     iterNT=',iterNT
        write(53,*)'Energy ',energy,' at iter ',iter,
     1    'kete(',jnr,')     iterNT=',iterNT
c note: all without dih4  of first CH3 group 
c   and     without dih22 of last  CH3 group
        rewind 14
        Read(14,*) (grad(i),i=1,nn)
        rewind 15
        Do 15  k=1,nn
15      Read(15,*) (H(k,i),i=1,nn)
        write(44,*) ' Diag of Hessian  '
        write(44,2) (H(I,I),I=1,5)
        write(44,2) (H(I,I),I=6,nn)
       endif
c c ccccccccccccccccccccccccccccccccccccccccccccccccc
c Ginv is the G-Matrix of Wilson !!! Ginv=B*B^T, g^{ij}
c Gmat=Ginv^{-1} is the usual covarinat metric matrix, g_{ij}
c all without mass-weighting
       rewind 48
       Do 16, k=1,nn
16       READ(48,*) (ginv(k,I),I=1,nn)
       READ(48,*)
       Do 17, k=1,nn
17       READ(48,*) (gmat(k,I),I=1,nn)
       READ(48,*)
c        IF(jnr.eq.1) then
c         write(44,*) ' ginv is from the B-mat '
c         Do 7711, k=1,nn
c7711     WRITE(44,2) (ginv(k,I),I=1,nn)
c         write(44,*) ' gmat '
c         Do 7712, k=1,nn
c7712     WRITE(44,2) (gmat(k,I),I=1,nn)
c        Endif
c testprint:  is o.k.
c        CALL matmult(ginv,nn,gmat,nn,unit,nn)
c        write(*,*)'test metric ginv*gmat for k=1 '
ccc        do k=1,nn
c             k=1
c             write(*,*)(unit(k,i),i=1,nn)
ccc        enddo
c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc       
c                                 use Proj for a help matrix
       call eigen(H,nn,nn,eival,eivec,Proj,0)
       write(44,*)'eigenvalues of  H '
       write(44,*)(eival(i),i=1,nn)
c       write(44,*)'eigenvectors of H leave out '
c       do k=1,nn
c          write(44,*)(eivec(k,i),i=1,nn)
c       enddo
       write(44,*)'first 10 eigenvector componentes phi psi  of H '
       do k=1,10
          write(44,*)k,'   ', eivec(k,20), eivec(k,44)
       enddo
       detHe=det(H,nn)
       write(44,*) 'det of  Hessian', detHe
       write(*,*)  'det of  Hessian', detHe
c c ccccccccccccccccccccccccccccccccccccccccccccccccc
        call vec_dinit(gradcont,nn,0.0d0)
        call matmult(ginv,nn,grad,1,gradcont,nn)
         write(44,*) ' grad '
         WRITE(44,2) (grad(I),I=1,5)
         WRITE(44,2) (grad(I),I=6,nn)
         write(44,*) ' gradcont '
         WRITE(44,2) (gradcont(I),I=1,5)
         WRITE(44,2) (gradcont(I),I=6,nn)
        rgnorm=dsqrt(vec_ddot(grad,gradcont,nn))
        WRITE(44,*) ' full gradient norm ',rgnorm
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Start uses coordinates-version of rsel from JTCC 09
c  rc for dyadic product is then covariant
c      if(iter.eq.0) then
c           call matmult(gmat,nn,rsel,1,rc,nn)
cc  or           
cc  turn grad-rsel (vri-dir) in contravariat version rc
cc  from part 3, rsel is covariant, like the gradient 
c      else
        call matmult(ginv,nn,rsel,1,rc,nn)
c      endif
cc  new:  projNm1: n-1 linearly independent lines from Proj
        CALL projector(Proj,rsel,rc,nn)
        write(44,*) '  projector by dyadic product '
c   old:  call ortcompl(rsel,projNm1,unit,nn)
c         write(44,*) ' projNm1 '
c         Do 7713, k=1,nn-1
c7713     WRITE(44,2) (projNm1(k,I),I=1,nn)
c  use n-1 lines of the dyadic product for the (n-1) x n matrix
c  of rank n-1
        do 19 k=1,nn-1
          do i=1,nn
            projNm1(k,i)=Proj(k,i)
          enddo
19      continue
cccccccccccccccccccccccccccccccccccccccccccccccccccc                     
c  print   
      if(jnr.eq.1) then
        write(44,*) '  rsel '
        WRITE(44,2) (rsel(I),I=1,5)
        WRITE(44,2) (rsel(I),I=6,nn)
        call vec_dinit(sw,nn,0.0d0)
        call matmult(projNm1,nn-1,rsel,1,sw,nn)
        write(44,*) ' Proj * rsel '
        WRITE(44,*) (sw(I),I=1,nn)
        call vec_dinit(sw,nn,0.0d0)
        call matmult(projNm1,nn-1,grad,1,sw,nn)
        write(44,*) ' Proj * grad '
        WRITE(44,*) (sw(I),I=1,nn)
      endif
cccccccccc WQ cccccccccccccccccccccccccccccccccccccccccccccccccc                 
c  proj: 1 x covariant, 1 x contravariant tensor               c
c  'Reduction' of Hesse matrix: 2x covariant                   c   
c   also rgeq is again two-fold covariant                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        call matmult(projNm1,nn-1,H,nn,rgeq,nn)
        write(44,*) 'print  rgeq leave out '
c        Do 7714, k=1,nn-1
c7714       WRITE(44,2) (rgeq(k,I),I=1,nn)
        call vec_dinit(rgrad,nn-1,0.0d0)
        call matmult(projNm1,nn-1,grad,1,rgrad,nn)
        write(44,*) ' rgrad '
        WRITE(44,2) (rgrad(I),I=1,nn-1)
        rgnorm=norm(rgrad,nn-1)
        WRITE(44,*) ' PES reduced grad norm  ',rgnorm
c Tangent ccccccccccccccccccccccccccccccccccccccccccccc
         call gettang(rgeq,gmat,tang,nn)
         write(44,*) 'tang from rgeq is '
         WRITE(44,2) (tang(I),I=1,5)
         WRITE(44,2) (tang(I),I=6,nn)
c       ccccccccccccccccccccccccccccccccccccccccccccc
        call matmult(rgeq,nn-1,tang,1,sw,nn)
        write(44,*) ' test: rgeq*tang gives '
        do k=1,nn/3
           write(44,*) (sw(3*(k-1)+I),I=1,3)
        enddo
        tgt=vec_mdot(tang,gmat,told,nn)
        write(44,*)' tang * gmat * told     ',tgt
        call vec_dinit(tcov,nn,0.0d0)
        call matmult(gmat,nn,tang,1,tcov,nn)
        write(44,*) '  tcov is '
        write(44,2) (tcov(I),I=1,5)
        write(44,2) (tcov(I),I=6,nn)
cccccccccccccccccccccccccccccccccccccccccccccc
        call buildk(dk,rgeq,tcov,rgrad,nn)
cccccccccccccccccccccccccccccccccccccccccccccc
c        write(44,*) ' K mat is '
c        Do 771, k=1,nn
c771       WRITE(44,*) (dk(k,I),I=1,nn+1)
        detk=det(dk,nn)
        write(*,*)  ' detK is       ', detk
        write(44,*) ' detK is       ', detk
        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_mdot(tang,gmat,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)
               call buildk(dk,rgeq,tcov,rgrad,nn)
          endif
          tg=vec_ddot(t1,t1old,nn)
          if (tg.lt.0.0d0 .and. jnr.gt.1 ) then 
             write(*,*)  ' bif of branches'
             write(44,*) ' bif of branches'
c  old: (does not work here) 
c             call vec_dmult_constant(tang,nn,-1.0d0,tang)
c             call vec_dmult_constant(tcov,nn,-1.0d0,tcov)
c             call buildk(dk,rgeq,tcov,rgrad,nn)
           endif
           call linsolve(dk,dstep,nn)
           WRITE(44,*) 'dstep by linsolve '
c           WRITE(44,2)(dstep(i),i=1,5)
c           WRITE(44,2)(dstep(i),i=6,nn)
cc        
ccccccc 1.order goto:
333      continue
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccc Abbruch 
        IF(
     1     rgnorm.lt.eps              .or.  ! Corrector o.k.
c    2        jnr.gt.1                .or.  ! predictor only
     3     iterNT.ge.iterNTmax        ) THEN
            write(44,*)' Predictor step ----------------  pstl=',pstl
            write(*,*) ' Predictor step ----------------  pstl=',pstl
            write(44,*) 'Nr j in vri2 is =',jnr
            write(*,*) '  |rgnorm|           ', rgnorm
            write(44,*) ' |rgnorm|           ', rgnorm
            iterNT=1
            rewind 17
            write(17,*) iter,iterNT
            rewind 25
            write(25,2)(tang(i),i=1,nn)
            rewind 26
            write(26,2)(t1(i),i=1,nn)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
           if(jnr.lt.le+1) THEN
c        c                               do predictorstep
            point(1)=kette(jnr,1)/btoa
            point(2)=kette(jnr,2)/btoa
            point(3)=kette(jnr,3)*torad
            point(4)=kette(jnr,4)/btoa
            point(5)=kette(jnr,5)*torad
          DO 20 k=2,18
            point(3*k  )= kette(jnr,3*k  )/btoa
            point(3*k+1)= kette(jnr,3*k+1)*torad
            point(3*k+2)= kette(jnr,3*k+2)*torad
20        continue
            point(57)=kette(jnr,57)/btoa
            point(58)=kette(jnr,58)*torad
cc predictor step by RGF 3
c            dk(nn,nn+1)=pstl
c            call linsolve(dk,dstep,nn)
c            call vec_dmult_constant(dstep,nn,0.25d0,dstep)
            call vec_dmult_constant(tang,nn,pstl,dstep)
            tgt=dmnorm(dstep,gmat,nn)
            write(*,*)
            write(44,*)'pred-step length ',tgt
            write(44,*)'pred-step by tangent would be'
            write(44,2)(dstep(i),i=1,5)
            write(44,2)(dstep(i),i=6,nn)
C C C C c C C C C c C C C C c C C C C
cc predictor step by Branin-Differential EQ
cc shift H to make a calculable determinant 
         write(44,*)'Hessian entries shifted by factor 4 for Ag'
         call vec_dmult_constant(H,nn*nn,4.0d0,H)
         detHe=det(H,nn)
         write(44,*) 'det of SHIFTED Hessian', detHe
ccc now: Adj calculable by the value-shift of H          
          call vec_dinit(Adj,nn*nn,0.0d0)  
          CALL mat_adj(H,nn,Adj)
          CALL matmult(Adj,nn,grad,1,Ag,nn)
          write(44,*)'Ag is  ...  '
          write(44,2)(Ag(i),i=1,5)
          write(44,2)(Ag(i),i=6,nn)
ccc turn Ag in direction of tang
          tgt=vec_mdot(tang,gmat,Ag,nn)
          write(44,*)' tang*gmat*Ag    ',tgt
          sig=1.0d0
          if (tgt.lt.0.0d0) sig=-1.0d0
          tminNT=dmnorm(Ag,gmat,nn) 
          factor=pstl/tminNT *sig 
          call vec_dmult_constant(Ag,nn,factor,dstep) 
          write(44,*)'pred-step length  set  to     ',pstl 
          write(44,*)'pred-step by Branin used '
          write(44,2)(dstep(i),i=1,5)
          write(44,2)(dstep(i),i=6,nn)
          write(*,*)'along NT |Ag|              = ',tminNT
     1                ,' at node ', jnr
          write(44,*)'along NT |Ag|             = ',tminNT
          write(44,*)'  at node ', jnr
c C  C C C C c c c c c c c c C C C C C c C C C C c C C C C c C C C C
             do i=1,nn
              point(i)=point(i)+dstep(i)
             enddo
c
             kette(jnr+1,1)=point(1)*btoa
             kette(jnr+1,2)=point(2)*btoa
             kette(jnr+1,3)=point(3)/torad
             kette(jnr+1,4)=point(4)*btoa
             kette(jnr+1,5)=point(5)/torad
            DO 21 k=2,18
             kette(jnr+1,3*k  )= point(3*k  )*btoa
             kette(jnr+1,3*k+1)= point(3*k+1)/torad
             kette(jnr+1,3*k+2)= point(3*k+2)/torad
21      continue
             kette(jnr+1,57)=point(57)*btoa
             kette(jnr+1,58)=point(58)/torad
c                     
       WRITE(44,*) 'ketten-point nach predictor phi & psi is '
       WRITE(44,*) '{', kette(jnr+1,20), ',', kette(jnr+1,44) ,' },'
       WRITE(7,*)  '{', kette(jnr+1,20), ',', kette(jnr+1,44) ,' },'
cccccccccccccccccccccccccccccccccccccccccccccccccc                                                                                
            jnr=jnr+1
            GOTO 55
c noch Abbruch                
       else
c           goto vrisearch in part Q3, or stop in Q5
                istatus=0
                rewind 18
                write(18,*) istatus
                icalc=0
                rewind 19
                write(19,*) icalc,sgn
                if(iter.eq.0) then
                   iter=1
                   rewind 17
                   write(17,*) iter,iterNT
c   drehe kette herum  
                   DO 23 k=1,le+1
                     do i=1,nn
                       rkette(k,i)=kette(le+2-k,i)
                     enddo
23                 continue     
                   DO 24 k=1,le+1
                     do i=1,nn
                       kette(k,i)=rkette(k,i)
                     enddo
24                 continue                          
                   rewind 25
                   write(25,2)(-1.0d0*tang(i),i=1,nn)
                   rewind 26
                   write(26,2)(-1.0d0*t1(i),i=1,nn)
                endif 
                goto 66
            ENDIF
ccc AbbruchEnde 
        ENDIF
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccCCCCCCCCCCCCCCCCCC
        tt=norm(dstep,nn)
        write(44,*) ' dstep norm Chartesian whould be ',tt
        tt=dmnorm(dstep,gmat,nn)
        write(44,*) ' dstep norm by gmat whould be    ',tt
        write(44,*) ' dstep is '
ccccc Notbremse  --  emergency break : 
      if(tt.gt.pstl)call vec_dmult_constant(dstep,nn,pstl/tt/2.,dstep)
        write(44,2) (dstep(I),I=1,5)
        write(44,2) (dstep(I),I=6,nn)
        write(44,*) ' dstep phi & psi is'
        write(44,2) dstep(20),dstep(44)
ccccc        
         DO 25 i=1,nn
             point(i)=point(i)+dstep(i)
25      CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
            WRITE(44,*) 'point nach dstep addition in Bohr/radian'
            WRITE(44,2) (point(I),I=1,5)
            WRITE(44,2) (point(I),I=6,nn)
c       
          kette(jnr,1)=point(1)*btoa
          kette(jnr,2)=point(2)*btoa
          kette(jnr,3)=point(3)/torad
          kette(jnr,4)=point(4)*btoa
          kette(jnr,5)=point(5)/torad
        DO 26 k=2,18
          kette(jnr,3*k  )= point(3*k  )*btoa
          kette(jnr,3*k+1)= point(3*k+1)/torad
          kette(jnr,3*k+2)= point(3*k+2)/torad
26      continue
          kette(jnr,57)=point(57)*btoa
          kette(jnr,58)=point(58)/torad
c          
          WRITE(44,*) 'ketten-point nach dstep addition phi & psi '
          WRITE(44,*) '{',kette(jnr,20),',',kette(jnr,44),' },' 
          WRITE(7,*)  '{',kette(jnr,20),',',kette(jnr,44),' },'
         write(44,*) 'Nr j in vri2 is =',jnr
         write(*,*)  'j vri2 = ',jnr,' iterNT  =  ',iterNT
         write(*,*) ' hypo |dstep|            ', tt,'   rg ',rgnorm
         write(44,*) 'hypo |dstep|            ', tt,'   rg ',rgnorm
               iterNT=iterNT+1
               rewind 17
               write(17,*) iter,iterNT               
55      continue
        rewind 9
        WRITE(9,*) jnr
66      continue
        if(icalc.eq.1) then
            WRITE(44,*) jnr
            WRITE(44,2) (kette(jnr,i),i=1,5) 
            WRITE(44,2) (kette(jnr,i),i=6,nn)
        endif
        rewind 8
        DO 29 k=1,le+1
29        WRITE(8,2)   (kette(k,i),i=1,nn)
c             
       rewind 18
       write(18,*) istatus
       rewind 19
       write(19,*) icalc,sgn
       call exit(istatus) 
      END
!
! Legt gerade "kette" der Laenge "le" zwischen die Punkte "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
cccccccccccccccccccccccccccccccccccccccccccccccc
! Berechnet den Projektionsoperator "pr" 
! aus der Suchrichtung  "r" covariant, und "rc" contravariat mit
!     pr= (Id-r*rc^T)  , pr is a mixed tensor      
!
       SUBROUTINE projector(pr,r,rc,nn)
        DOUBLE PRECISION r(nn),rc(nn),pr(nn,nn)
        DOUBLE PRECISION rn(nn),rcn(nn),nrm,vec_ddot
c        nrm=1.d0*norm(r,nn)
        nrm=dsqrt(vec_ddot(r,rc,nn))
        DO 5 i=1,nn
         rn(i)=r(i)/nrm
         rcn(i)=rc(i)/nrm
5     CONTINUE
      DO 20 i=1,nn
       DO 10 j=1,nn
        IF (i.EQ.j) THEN
          pr(i,j)=-rn(i)*rcn(j)+1.d0
        ELSE
          pr(i,j)=-rn(i)*rcn(j)
        ENDIF
10      CONTINUE
20    CONTINUE
      END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccc                                       cccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine ginvers(in,NAT,out,omat_inv)
c    calculates the inverse matrix of IN
      real*8 in,out
      logical omat_inv
      integer NAT,i,j,k
      dimension in(NAT,NAT),out(NAT,NAT)
      do 1,k=1,NAT
      do 1,j=1,NAT
  1   out(j,k)=in(j,k)
c
      do 10,k=1,NAT
         if (out(k,k).eq.0.0d0) then
            omat_inv=.false.
            return
         endif
         out(k,k)=1.0d0/out(k,k)
         do 20,i=1,NAT
             if (i.ne.k) out(i,k)=-(out(i,k)*out(k,k))
   20    continue
         do 30,j=1,NAT
             if (j.ne.k) then
                do 40,i=1,NAT
                      if (i.ne.k) out(i,j)=out(i,j)+out(i,k)*out(k,j)
   40           continue
                out(k,j)=out(k,j)*out(k,k)
             endif
   30    continue
   10 continue
      omat_inv=.true.
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       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
cccccccccccccccccccccccccccccccc
!
! Multiplikation Ax=b
!
      SUBROUTINE multiply(b,A,x,nn)
      DOUBLE PRECISION b(nn),A(nn,nn),x(nn)
      DO 20 i=1,nn
        b(i)=0
        DO 10 j=1,nn        
           b(i)=b(i)+A(i,j)*x(j)
10      CONTINUE
20    CONTINUE
      END
!
! 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-23) THEN
        norm=0.d0
      ELSE
        norm=DSQRT(s)
      ENDIF
      END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
       subroutine ortcompl(vec,compl,dmet,ndim)
c   calculate an Orthonormalcomplement COMPL to vector VEC
       real*8 basis,vec,compl,swap,det,dmgs,dmet
       integer i,j,ndim
       dimension vec(ndim),compl(ndim-1,ndim),basis(ndim,ndim),
     *   swap(ndim),dmet(ndim,ndim)
       call vec_dinit(basis,ndim*ndim,0.0d0)
       do 10,i=1,ndim
         basis(i,i)=1.0d0
 10    continue
       do 20,i=1,ndim
         if (vec(i).ne.0.0d0) then
           do 30,j=1,ndim
             swap(j)=basis(j,1)
             basis(j,1)=vec(j)
             if (i.ne.1) basis(j,i)=swap(j)
 30        continue
           goto 40
         endif
 20    continue
 40    continue
c  note: dmgs is the Gram-Schmidt method 
       det=dmgs(basis,dmet,ndim)
       call vec_dinit(compl,(ndim-1)*ndim,0.0d0)
       call mat_trans(basis,ndim,ndim,basis)
       do 60,i=2,ndim
          call vec_mrowcopy(i,basis,ndim,i-1,compl,ndim-1,ndim)
 60    continue
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       real*8 function dmgs(a,dmet,dim)
c  modified Gram-Schmidt algorithm return det(a)
       real*8 det,rr,a,q,r,vec_mmdot,dmet,dzero
       integer dim,i,k
       dimension a(dim,dim),q(dim,dim),r(dim,dim),dmet(dim,dim)
       data dzero/0.0d0/
       call vec_dinit(q,dim*dim,0.0d0)
       call vec_dinit(r,dim*dim,0.0d0)
       det=1.0d0
       call mat_dcopy(a,q,dim,dim)
       do 10,k=1,dim
         if (k.eq.1) goto 21
         do 20,i=1,k-1
            r(i,k)=vec_mmdot(i,q,dim,k,q,dim,dim,dmet)
            call vec_madd(k,q,dim,1.0d0,i,q,dim,-r(i,k),k,q,dim,dim)
 20      continue
 21      continue
         rr=vec_mmdot(k,q,dim,k,q,dim,dim,dmet)
         r(k,k)=dsqrt(rr)
         det=det*r(k,k)
         if(rr.lt.dzero) goto 33
         call vec_madd(k,q,dim,1.0d0/r(k,k),k,q,dim,0.0d0,k,q,dim,dim)
 10    continue
       call mat_dcopy(q,a,dim,dim)
       dmgs=det
       return
 33    CONTINUE
cc       call cerr(-23)
       return
       end
c
      subroutine gettang(curveq,gmat,tangent,nn)
c   This routine calculate the tangent on a curve
c   whereas the tangent is unique determined by:
c
c   H*t=0
c   |t|=1
c   det[H,t]>0
c
      real*8 tangent,gmat,curveq,q,ss1,det,zero,dmnorm
      integer nn,i,j
      dimension tangent(nn),curveq(nn-1,nn),q(nn*nn),ss1(nn,nn)
      dimension gmat(nn,nn)
      data zero /0.0d0/
      call qr(curveq,q,nn-1)
      call vec_dcopy(q(nn*(nn-1)+1),tangent(1),nn)
c test     do 1,i=1,nn
c  1   tangent(i)=q(nn,i)
c wq
       call vec_dmult_constant
     +   (tangent,nn,1.0d0/dmnorm(tangent,gmat,nn),tangent)
cc       call mat_dcopy(curveq,ss1,nn-1,nn)
       do 10,j = 1,nn-1
       do 20,i = 1,nn
          ss1(j,i) = curveq(j,i)
  20    continue
  10    continue
       do 30,i = 1,nn
  30      ss1(nn,i)=tangent(i)
       if (det(ss1,nn).lt.zero) then
        call vec_dmult_constant(tangent,nn,-1.0d0,tangent)
      endif
      return
      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,xx
      integer nn
      dimension vec(nn),gm(nn,nn)
      xx= vec_mdot(vec,gm,vec,nn)
c      if(dabs(xx) .lt. 1.0d-63) then
c              xx= 1.0d-63
c      endif
      dmnorm=dsqrt(dabs(xx))
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine qr(gex,q,n)
C      QR-Zerlegung von gex
       integer i,j,k,n
       real*8 s1,s2,s,gex,q,r,qq,rr
       dimension gex(n,n+1),q(n+1,n+1),r(n+1,n),
     * qq(n+1,n+1),rr(n+1,n)
       do 10,i=1,n
          do 20,j=1,n+1
             r(j,i)=gex(i,j)
 20       continue
 10    continue
       do 30,i=1,n+1
          do 40,j=1,n+1
             q(i,j)=0.0d0
 40       continue
          q(i,i)=1.0d0
 30    continue
       do 50,i=1,n
          do 60,j=i+1,n+1
             s1=r(i,i)
             s2=r(j,i)
             s=dsqrt(s1*s1+s2*s2)
             if (s.gt.0.0d0) then
                s1=s1/s
                s2=s2/s
                else
                s1=0.0d0
                s2=0.0d0
                endif
             call vec_dcopy(r,rr,(n+1)*n)
             do 70,k=1,n
                rr(i,k)=s1*r(i,k)+s2*r(j,k)
                rr(j,k)=s1*r(j,k)-s2*r(i,k)
 70          continue
             call vec_dcopy(rr,r,(n+1)*n)
             call vec_dcopy(q,qq,(n+1)*(n+1))
             do 80,k=1,n+1
                qq(i,k)=s1*q(i,k)+s2*q(j,k)
                qq(j,k)=s1*q(j,k)-s2*q(i,k)
 80          continue
             call vec_dcopy(qq,q,(n+1)*(n+1))
 60       continue
 50    continue
       do 90,i=1,n+1
          do 100,j=1,n+1
             q(i,j)=qq(j,i)
 100      continue
 90    continue
c       write(44,*) 'in QR '
c       do i=1,n+1
c         write(44,*) (q(i,j),j=1,n+1)
c       enddo
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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 buildk(dk,req,tang,rgr,nn)
c build the K-Matrix
      real*8 dk,req,tang,rgr
      integer nn,i,j
      dimension dk(nn,nn+1),req(nn-1,nn),tang(nn),rgr(nn-1)
      do 10,i=1,nn-1
        do 11,j=1,nn
        dk(i,j)=req(i,j)
 11   continue
 10   continue
      do 20,j=1,nn
         dk(nn,j)=tang(j)
 20   continue
      do 30,i=1,nn-1
         dk(i,nn+1)=-rgr(i)
 30   continue
      dk(nn,nn+1)=0.0d0
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
       subroutine linsolve(matein,x,n)
c  solves MATEIN * X = 0
       real*8 matein,dre,x,w
       integer n
       dimension matein(n,n+1),dre(n,n+1),x(n),w(n)
       call householder(matein,dre,w,n,n+1)
c       write(44,*) 'in linsolve '
c       do k=1,n
c         write(44,*) (dre(k,i),i=1,n+1)
c       enddo
       call loesung(dre,x,n)
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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 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
c
       subroutine loesung(matein,x,n)
c   intern linear equation routine
       real*8 matein,x
       integer i,j,n,k
       dimension matein(n,n+1),x(n)
       do 10,i=1,n
 10     x(i)=matein(i,n+1)
       x(n)=x(n)/matein(n,n)
       do 20,i=1,n-1
         k=n-i
         do 30,j=k+1,n
 30        x(k)=x(k)-matein(k,j)*x(j)
           x(k)=x(k)/matein(k,k)
 20    continue
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
       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
ccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine vec_dmult_add(arra1,arra2,ndim,zahl,arra4)
c  For each Komponent: ARRA4 = ARRA1 + ZAHL*ARRA2
      real*8  arra1,arra2,arra4,zahl
      integer ndim,i
      dimension arra1(ndim),arra2(ndim),arra4(ndim)
        if(ndim .le. 0)return
           do 145 i = 1,ndim
             arra4(i) = arra1(i)+arra2(i)*zahl
 145       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
      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
       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
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine eigen (a,m,n,d,vec,e,iff)
c  eigenvalues and -vectors of A      
      real*8 a,d,vec,e,eps,tol,h,g,s,f,b,p,r,c
      integer m,iff,n,i,j,ni,ii,l,k,j1
      dimension a(m,m), d(m), vec(m,m)
      dimension e(m)
      eps = 0.278d-16
      tol = 0.211758d-21
      if (n.eq.1) then
         d(1) = a(1,1)
         vec(1,1) = 1.0d0
      else
         do 30 i = 1 , n
c     householder's reduction
c     simulation of loop do 150 i=n,2,(-1)
            do 20 j = 1 , i
               vec(i,j) = a(i,j)
 20         continue
 30      continue
         do 120 ni = 2 , n
            ii = n + 2 - ni
            do 110 i = ii , ii
               l = i - 2
               h = 0.0d0
               g = vec(i,i-1)
               if (l.gt.0) then
                  do 40 k = 1 , l
                     h = h + vec(i,k)**2
 40               continue
                  s = h + g*g
                  if (s.lt.tol) then
                     h = 0.0d0
                  else if (h.gt.0) then
                     l = l + 1
                     f = g
                     g = dsqrt(s)
                     if (f.gt.0) then
                        g = -g
                     end if
                     h = s - f*g
                     vec(i,i-1) = f - g
                     f = 0.0d0
                     do 70 j = 1 , l
                        vec(j,i) = vec(i,j)/h
                        s = 0.0d0
                        do 50 k = 1 , j
                           s = s + vec(j,k)*vec(i,k)
 50                     continue
                        j1 = j + 1
                        if (j1.le.l) then
                           do 60 k = j1 , l
                              s = s + vec(k,j)*vec(i,k)
 60                        continue
                        end if
                        e(j) = s/h
                        f = f + s*vec(j,i)
 70                  continue
                     f = f/(h+h)
                     do 80 j = 1 , l
                        e(j) = e(j) - f*vec(i,j)
 80                  continue
                     do 100 j = 1 , l
                        f = vec(i,j)
                        s = e(j)
                        do 90 k = 1 , j
                           vec(j,k) = vec(j,k) - f*e(k) - vec(i,k)*s
 90                     continue
 100                 continue
                  end if
               end if
c     accumulation of transformation matrices
               d(i) = h
               e(i-1) = g
 110        continue
 120     continue
         d(1) = vec(1,1)
         vec(1,1) = 1.0d0
         do 170 i = 2 , n
            l = i - 1
            if (d(i).gt.0) then
               do 150 j = 1 , l
                  s = 0.0d0
                  do 130 k = 1 , l
                     s = s + vec(i,k)*vec(k,j)
 130              continue
                  do 140 k = 1 , l
                     vec(k,j) = vec(k,j) - s*vec(k,i)
 140              continue
 150           continue
            end if
            d(i) = vec(i,i)
            vec(i,i) = 1.0d0
            do 160 j = 1 , l
c     diagonalization of the tridiagonal matrix
               vec(i,j) = 0.0d0
               vec(j,i) = 0.0d0
 160        continue
 170     continue
         b = 0.0d0
         f = 0.0d0
         e(n) = 0.0d0
         do 250 l = 1 , n
c     test for splitting
            h = eps*(dabs(d(l))+dabs(e(l)))
            if (h.gt.b) b = h
            do 180 j = l , n
c     test for convergence
               if (dabs(e(j)).le.b) go to 190
 180        continue
 190        if (j.eq.l) then
               d(l) = d(l) + f
            else
 200           p = (d(l+1)-d(l))*0.5d0/e(l)
               r = dsqrt(p*p+1.0d0)
               if (p.lt.0) then
                  p = p - r
               else
                  p = p + r
               end if
               h = d(l) - e(l)/p
               do 210 i = l , n
c     qr transformation
                  d(i) = d(i) - h
 210           continue
               f = f + h
               p = d(j)
c     simulation of loop do 330 i=j-1,l,(-1)
               c = 1.0d0
               s = 0.0d0
               j1 = j - 1
               do 240 ni = l , j1
                  ii = l + j1 - ni
                  do 230 i = ii , ii
c     protection against underflow of exponents
                     g = c*e(i)
                     h = c*p
                     if (dabs(p).lt.dabs(e(i))) then
                        c = p/e(i)
                        r = dsqrt(c*c+1.0d0)
                        e(i+1) = s*e(i)*r
                        s = 1.0d0/r
                        c = c/r
                     else
                        c = e(i)/p
                        r = dsqrt(c*c+1.0d0)
                        e(i+1) = s*p*r
                        s = c/r
                        c = 1.0d0/r
                     end if
                     p = c*d(i) - s*g
                     d(i+1) = h + s*(c*g+s*d(i))
                     do 220 k = 1 , n
                        h = vec(k,i+1)
                        vec(k,i+1) = vec(k,i)*s + h*c
                        vec(k,i) = vec(k,i)*c - h*s
 220                 continue
 230              continue
 240           continue
               e(l) = s*p
c     convergence
               d(l) = c*p
               if (dabs(e(l)).gt.b) go to 200
c     ordering of eigenvalues
               d(l) = d(l) + f
            end if
 250     continue
         if (iff.lt.1) then
            ni = n - 1
            do 280 i = 1 , ni
               k = i
               p = d(i)
               j1 = i + 1
               do 260 j = j1 , n
                  if (d(j).lt.p) then
                     k = j
                     p = d(j)
                  end if
 260           continue
               if (k.ne.i) then
                  d(k) = d(i)
                  d(i) = p
                  do 270 j = 1 , n
                     p = vec(j,i)
                     vec(j,i) = vec(j,k)
                     vec(j,k) = p
 270              continue
               end if
c     special treatment of case n = 1
 280        continue
         end if
      end if
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
      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),
     1          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-70) goto 20
                matout(i,j)= (-1.0d0)**(i+j)*matout(i,j)
 20     continue
 10   continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccc
