      PROGRAM HCNvri2
!
!   par t(2) rgf to dir=rsel  
! 
!  W.Quapp, Benjamin Schmidt, Decembre 2009
!
!   nn     = Dimension,  N3 = number of atoms
!   le     = Chainlength (even integer)
! iterNTmax= maximal number of NT correctors: if corector
!            uses more than, say 4, steps: it works not fine.. 
! pst,pstl = Steplength of predictor, depends from iteration 
!    eps   = Convergence in RGF corrector, red from part 1
!                      is useless for predictors only
!  kette   = Chain between points xstart and xend,
!                      xend is useless for predictors only
!   rsel   = 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 part:
! (2) To "dir" calculate corresponding NT, by pred-corr,
!              or by predictors along AG only
!
!     lather parts:
! (3) Calculate VRI approximation,
!     one version: use its gradient for new "dir"
!     other version: use only Ag for a predictor step
! (4) Organize points for Gamess-US input for (2) or (3)
! HCNvriReadOUT: 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 stl,pstl
      PARAMETER (nn=3,N3=3,le=50,iterNTmax=10,stl=1.d0,pst=0.012d0)
!
! RGF-Verfahren
!
      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,
     7  t1old(nn),t1(nn),rgnorm,point(nn),eival(nn),eivec(nn,nn),
     8  gradcont(nn),ginv(nn,nn),energy,sgn
      INTEGER k,j,i,jnr,istatus,icalc,iter,iterNT,Natoms(N3)
      logical omat_inv
      Character*2  CharAtoms(N3),Zcoor(nn)
      Character*70 zmatLines(20)
       Data Zcoor/'hc','nc','w '/
       Data CharAtoms/'H ','C ','N '/
       Data Natoms/1,6,7/
c     pure and applied chemistry, 51, 1 (1979):
c     angstroms per bohr
       data toang  /0.52917706d0/
ccccccccccccccccccccccccccccccccccccccccc      
      OPEN(8,FILE='hcnVkette.dat')
c      OPEN(91,FILE='hcnVkette.mol')
c      OPEN(95,FILE='hcnVkette.gau')
ccccccccccccccccccccccccccccccccccccccccc
      OPEN(7,FILE='chain.dat',status='unknown',access='append')
      OPEN(9,FILE='hcnVketNrj.dat')
      OPEN(11,FILE='hcnVdir.dat')
      OPEN(13,FILE='hcnVener.dat')
      OPEN(53,FILE='hcnVEall.dat',status='unknown',access='append')
      OPEN(14,FILE='hcnVgrad.dat')
      OPEN(15,FILE='hcnVhesse.dat')
      OPEN(16,FILE='hcnVeps.dat')
      OPEN(17,FILE='hcnViter.dat')
      OPEN(18,FILE='hcnVistat.dat')
      OPEN(19,FILE='hcnVicalc.dat')
      OPEN(20,FILE='hcnVketst.dat')
      OPEN(21,FILE='hcnVketen.dat')
      OPEN(25,FILE='hcnVtold.dat')
      OPEN(26,FILE='hcnVt1old.dat')
      OPEN(28,FILE='hcnVketNrj1.dat')
      OPEN(32,FILE='hcnVnode.dat')
      OPEN(33,FILE='hcnVmin.dat')
      OPEN(48,FILE='hcnVgmatr.dat')
c protocol file: all interesting output of the vri-search
      OPEN(44,FILE='hcnVprotocol.txt',status='unknown',access='append')
cccccccccccccccccccccccccccccccccccccccc
2     FORMAT(3F20.14)
      N=nn
c     j1=le/2
      j1=1
      k1=1  
c     k1=le+1
c     j1,k1 for iteration in part 3
      rewind 28
      write(28,*) j1,k1 
      rewind 16
      read(16,*) eps
c  eps is criterion for rgf convergence to given rsel-NT
       rewind 9
       Read(9,*) jnr
c jnr is the currend node number
      If(jnr.eq.2) then
         write(*,*) 'eps = ',eps
         write(44,*)'eps = ',eps
       endif
 
       rewind 8
       DO 1 k=1,le+1
          Read(8,*) (kette(k,i),i=1,nn)
1      CONTINUE
 
      rewind 19
      read(19,*) icalc,sgn
      rewind 17    
      Read(17,*) iter,iterNT
      if(iter.eq.0) then
        pstl= pst
        else
        pstl= 2.0d0*pst
      endif   
      if(iterNT.gt.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,') in angstroem,degree HCNvri2 '
         Do 7, i=1,nn 
7         write(44,*) Zcoor(i), 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
c Note:
c in corrector loop with use of gradient and hessian, there
c are used bohr and radian. Input of Gamess, however, is in
c angstroem and degree !!!  
c c ccccccccccccccccccccccccccccccccccccccccccccccccc
      call mat_diag(unit,nn,1.0d0)
      rewind 11
      Read(11,*) (rsel(i),i=1,nn)
      if(jnr.eq.2) then  
        write(*,*)  'dir in vri2  ', (rsel(i),i=1,nn)
        write(44,*) 'dir in vri2  ', (rsel(i),i=1,nn)
      endif
      rewind 18 
      read(18,*) istatus
c
      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              
       rewind 19
       read(19,*) icalc,sgn
       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
        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)
       endif
cc           use  'Proj' for a help matrix       
c       call eigen(H,nn,nn,eival,eivec,Proj,0)
c       write(44,*)'eigenvalues of H '
c       write(44,*)(eival(i),i=1,nn)
c       write(44,*)'eigenvectors of H '
c       do k=1,nn
c          write(44,*)(eivec(k,i),i=1,nn)
c       enddo
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 note: all is without mass-weighting
        rewind 48
        Do 766, k=1,nn
766     READ(48,*) (ginv(k,I),I=1,nn)
        READ(48,*)
        Do 767, k=1,nn
767     READ(48,*) (gmat(k,I),I=1,nn)
        READ(48,*)
        IF(jnr.eq.1) then
         write(44,*) ' ginv is from the B-mat '
         Do 7711, k=1,nn
7711     WRITE(44,2) (ginv(k,I),I=1,nn)
         write(44,*) ' gmat '
         Do 7712, k=1,nn
7712     WRITE(44,2) (gmat(k,I),I=1,nn)
        Endif
c testprint
c        CALL matmult(ginv,nn,gmat,nn,unit,nn)
c            write(*,*)'test metric ginv*gmat'
c            do k=1,nn
c              write(*,*)(unit(k,i),i=1,nn)
c            enddo
c unit, gmat, ginv is for Metric in internal coordinates 
        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,nn)
         write(44,*) ' gradcont '
         WRITE(44,2) (gradcont(I),I=1,nn)
        gnorm=dsqrt(vec_ddot(grad,gradcont,nn))
        WRITE(44,*) ' PES full gradient norm ',gnorm
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  rsel has to be covariant, like the gradient 
c  projNm1: n-1 linearly independent lines from Proj  
      call matmult(ginv,nn,rsel,1,rc,nn)
c     write(44,*) 'dir contravariant ',(rc(i),i=1,nn)
      CALL projector(Proj,rsel,rc,nn)
c old, is ruled out:  call ortcompl(rsel,projNm1,unit,nn)
c        write(44,*) ' Proj '
c        Do 7713, k=1,nn
c7713     WRITE(44,2) (Proj(k,I),I=1,nn)
        do 7719 k=1,nn-1
          do i=1,nn
            projNm1(k,i)=Proj(k,i)
          enddo
7719    continue
cccccccccccccccccccccccccccccccccccccccccccccccccccc
c test print :  is o.k. 
c        write(44,*) '  rsel '
c        WRITE(44,*) (rsel(I),I=1,nn)
c        call matmult(proj,nn,rsel,1,sw,nn)
c        write(44,*) ' projNm1 * rsel '
c        WRITE(44,*) (sw(I),I=1,nn)
cccccccccc WQ   2009 cccccccccccccccccccccccccccccccccccccc
c                                                         c
c  proj: 1 x covariant, 1 x contravariant tensor          c
c  'Reducting' of Hessian:  is again 2x covariant         c
c        rgeq is again two-fold covariant                 c
c                                                         c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        call matmult(projNm1,nn-1,H,nn,rgeq,nn)
        write(44,*) ' rgeq '
        Do 7714, k=1,nn-1
7714       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)
c sw: use full dimension of rgrad for metric
        call matmult(Proj,nn,grad,1,sw,nn)
        write(44,*) ' rgrad '
        WRITE(44,*) (sw(I),I=1,nn)
        rgnorm=dmnorm(sw,ginv,nn)
        WRITE(44,*) ' PES reduced grad norm  ',rgnorm
c Tangent ccccccccccccccccccccccccccccccccccccccccccccc
           call gettang(rgeq,gmat,tang,nn)
           write(44,*) 'tang from rgeq used '
           WRITE(44,*) (tang(I),I=1,nn)
c near VRI it is better to use Ag for the tangent
c        CALL mat_adj(H,nn,Adj)
c        CALL matmult(Adj,nn,grad,1,sw,nn)
c        tgt=dmnorm(sw,gmat,nn)
c        do i=1,nn
c          sw(i)=sw(i)/(-1.0d0*tgt)
c        enddo  
c        write(44,*) 'tang by Ag'
c        WRITE(44,*) (sw(I),I=1,nn)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        call matmult(rgeq,nn-1,tang,1,sw,nn)
        write(44,*) 'rgeq*tang gives ',(sw(I),I=1,nn-1)
        tgt=vec_mdot(tang,gmat,told,nn)
        write(44,*)'tang*G_ij*told     ', tgt
        call vec_dinit(tcov,nn,0.0d0)
        call matmult(gmat,nn,tang,1,tcov,nn)
        write(44,*) '  tang is '
        write(44,*) (tang(I),I=1,nn)
        write(44,*) '  tcov is '
        write(44,2) (tcov(I),I=1,nn)
cccccccccccccccccccccccccccccccccccccccccccccc
        call buildk(dk,rgeq,tcov,rgrad,nn)
cccccccccccccccccccccccccccccccccccccccccccccc
        write(44,*) ' K mat is '
        Do 7715, k=1,nn
7715       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
cccccccccccccccccccccccccccccccccccccccccccccc
          tg=vec_ddot(t1,t1old,nn)
          if (tg.lt.0.0d0 .and. jnr.gt.1 ) then 
             write(*,*)  ' bif of branches'
             write(44,*) ' bif of branch
c   that is ruled out ??
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 would be '
           WRITE(44,*)(dstep(i),i=1,nn)
c
c  cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  ccccccc 1.order goto:
333      continue
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccc Abbruch 
c                 Abbruch --  condition start
c       IF(jnr.gt.4  ) then 
         write(44,*) ' Predictor steps only *************  pstl=',pstl
         write(*,*)  ' Predictor steps only *************  pstl=',pstl
c       endif
c       IF(
cc      1    Dabs(detk).lt.0.0005d0      .or. 
c     1     rgnorm.lt.eps              .or.  
c     2     jnr.gt.4                   .or.
c     3     iterNT.gt.10             ) THEN
ccc          write(44,*) ' Predictor steps only **************** '
ccc          write(*,*) ' Predictor steps only **************** '
            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
cc do predictorstep
            point(1)=kette(jnr,1)/btoa
            point(2)=kette(jnr,2)/btoa
            point(3)=kette(jnr,3)*torad
cc test predictor step by RGF 3
            dk(nn,nn+1)=pstl
            call linsolve(dk,dstep,nn)
            tgt=dmnorm(dstep,gmat,nn)
            write(*,*)
            write(44,*)'pred-step length',tgt
            write(44,*)'pred-step by RGF3 would be'
            write(44,*)(dstep(i),i=1,nn)
C C C C c C C C C c C C C C c C C C C
c predictor step by Branin differential equ
          rewind 33
          read(33,*) tmin, Agmin, AgmiNT
          CALL mat_adj(H,nn,Adj)
          CALL matmult(Adj,nn,grad,1,Ag,nn)
          tminNT=dmnorm(Ag,gmat,nn) 
          write(44,*) 'tminNT of Ag  ',tminNT
          factor=pstl/tminNT 
          tgt=vec_mdot(Ag,gmat,told,nn)
        IF(jnr.lt.30 ) then
          if(tgt.lt. 0.0d0) factor=-1.0d0*factor
        endif    
        IF(jnr.ge.30 )factor=-1.0d0*factor
        call vec_dmult_constant(Ag,nn,factor,dstep) 
        tgt=dmnorm(dstep,gmat,nn)
        write(44,*)'pred-step length',tgt
        write(44,*)'pred-step by Branin is'
        write(44,*)(dstep(i),i=1,nn)
c additional: search lowest node for |Ag|
        IF(tminNT .LT. AgmiNT ) THEN
          rewind 33
          write(33,*) tmin, Agmin, tminNT
        ENDIF
           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
             do i=1,nn
              point(i)=point(i)+dstep(i)
             enddo
             kette(jnr+1,1)=point(1)*btoa
             kette(jnr+1,2)=point(2)*btoa
             kette(jnr+1,3)=point(3)/torad
             rewind 32
             write(32,*) (kette(jnr+1,i),i=1,nn)
cccccccccccccccccccccccccccccccccccccccccccccccccc                                                                                
            jnr=jnr+1
                GOTO 55
            else
c goto vrisearch in part 3, or stop in part 5
                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
                   DO 891 k=1,le+1
                     do i=1,nn
                       rkette(k,i)=kette(le+2-k,i)
                     enddo
891                continue     
                   DO 892 k=1,le+1
                     do i=1,nn
                       kette(k,i)=rkette(k,i)
                     enddo
892                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 Abbruch --  condition end 
c        ENDIF
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccCCCCCCCCCCCCCCCCCC
        tt=dmnorm(dstep,gmat,nn)
        write(44,*) ' dstep is '
ccccc Notbremse  --  emergency break : 
        factor=1.0d0 
        if(tt .gt. 0.05d0) factor=0.05d0/tt
        write(44,2) (dstep(I)*factor*stl,I=1,nn)
        DO 20 i=1,nn
           point(i)=point(i)+dstep(i)*factor*stl
20      CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        WRITE(44,*) 'point nach dstep addition'
        WRITE(44,*) (point(I),I=1,nn)
        kette(jnr,1)=point(1)*btoa
        kette(jnr,2)=point(2)*btoa
        kette(jnr,3)=point(3)/torad
        write(44,*) 'Nr j in vri2 is =',jnr
        write(*,*)  'j vri2 = ',jnr,' iterNT  =  ',iterNT
        write(*,*) '  |dstep|                ', tt,'   rg ',rgnorm
        write(44,*) ' |dstep|                ', tt,'   rg ',rgnorm
        iterNT=iterNT+1
        rewind 17
        write(17,*) iter,iterNT               
ccc
55      continue
ccc
        rewind 9
        WRITE(9,*) jnr
ccc
66      continue
ccc
202     FORMAT(I5,3X,3F20.14)
        if(icalc.eq.1) then
        DO 88 k=1,le+1
88        WRITE(44,202) k, (kette(k,i),i=1,nn)
        endif
        rewind 8
        DO 89 k=1,le+1
89      WRITE(8,2)       (kette(k,i),i=1,nn)
        rewind 32
        WRITE(32,2)       (kette(jnr,i),i=1,nn)
       if(jnr.eq.le+1  .and. iterNT.eq.1  )  then
c make Mma input file !! CH and CN coords interchanged
          WRITE(7,*)'{   (*  iter = ',iter,'  *) '
          DO 606 k=1,le
606       WRITE(7,632) kette(k,2),kette(k,1),kette(k,3)
          k=le+1
          WRITE(7,633) kette(k,2),kette(k,1),kette(k,3)
          WRITE(7,*)'}'
632       FORMAT('{',F20.14,' , ',F20.14,'  , ',F20.14,'/90.0 },')
633       FORMAT('{',F20.14,' , ',F20.14,'  , ',F20.14,'/90.0 }' )
        endif
       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 nrm,vec_ddot
c        nrm=1.d0*norm(r,nn)
        nrm=dsqrt(vec_ddot(r,rc,nn))
        DO 5 i=1,nn
         r(i)=r(i)/nrm
         rc(i)=rc(i)/nrm
5     CONTINUE
      DO 20 i=1,nn
       DO 10 j=1,nn
        IF (i.EQ.j) THEN
          pr(i,j)=-r(i)*rc(j)+1.d0
        ELSE
          pr(i,j)=-r(i)*rc(j)
        ENDIF
10      CONTINUE
20    CONTINUE
      END
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
Ccccccccccccccccccccccccccccccccccccccccccccc
! Berechnet adjunkte Matrix A zu M for n=3
      SUBROUTINE adjoint(A,M,nn)
      DOUBLE PRECISION A(nn,nn),M(nn,nn)
      A(1,1)=M(2,2)*M(3,3)-M(2,3)*M(3,2)
      A(1,2)=M(1,3)*M(3,2)-M(1,2)*M(3,3)
      A(1,3)=M(1,2)*M(2,3)-M(1,3)*M(2,2)
      A(2,1)=M(2,3)*M(3,1)-M(2,1)*M(3,3)
      A(2,2)=M(1,1)*M(3,3)-M(1,3)*M(3,1)
      A(2,3)=M(1,3)*M(2,1)-M(1,1)*M(2,3)
      A(3,1)=M(2,1)*M(3,2)-M(2,2)*M(3,1)
      A(3,2)=M(1,2)*M(3,1)-M(1,1)*M(3,2)
      A(3,3)=M(1,1)*M(2,2)-M(1,2)*M(2,1)
c       case n=2 
c       A(1,1)=M(2,2)
c       A(1,2)=-M(1,2)
c       A(2,1)=-M(2,1)
c       A(2,2)=M(1,1)
      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
ccccccccccccccccccccccccccccccccccccccc
! 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
       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
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine gettang(curveq,gmat,tangent,nn)
c   This routine calculate the tangent on a curve
c   whereas the tangent is unique determined by:
c   H*t=0
c   |t|=1
c   det[H,t]>0
      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)
       call vec_dmult_constant
     +   (tangent,nn,1.0d0/dmnorm(tangent,gmat,nn),tangent)
       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
      integer nn
      dimension vec(nn),gm(nn,nn)
      dmnorm=dsqrt(vec_mdot(vec,gm,vec,nn))
      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
       write(44,*) 'in QR '
       do i=1,n+1
         write(44,*) (q(i,j),j=1,n+1)
       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
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
       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
ccccccccccccccccccccccccccccccccccccccccccc
      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
cccccccccccccccccccccccccccccccccccccc
       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)
       write(44,*) 'in linsolve '
       do k=1,n
         write(44,*) (dre(k,i),i=1,n+1)
       enddo
       call loesung(dre,x,n)
       return
       end
cccccccccccccccccccccccccccccccccccccccccccc
      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
C            w(k)= 0.0d0
C            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
ccccccccccccccccccccccccccccccccccccc
       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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       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
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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       
ccccccccccccccccccccccccccccccccccccccccccccccc
       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
cccccccccccccccccccccccccccccccccccccccccccc
       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
cccccccccccccccccccccccccccccccccccccccccccccccccc
      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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       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-25) then
          det=0.0d0
          goto 20
       endif
       if (dabs(det) .lt.1.0d-25) then
          det=0.0d0
          goto 20
       endif
 10    det=det*mataus(i,i)
       if (iand(nn,1).eq.0) det=-det 
 20    continue
       return
       end
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
ccccccccccccccccccccccccccccccccccccccccccc
      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)
c                       if(dabs(matscr(ii,jj)).lt.1.0d-30) then
c                         matscr(ii,jj)=0.0d0
c                       endif
 40               continue
 30             continue
                matout(i,j)=det(matscr,ndim-1)
                if(dabs(matout(i,j)).lt.1.0d-30) goto 20
                matout(i,j)= (-1.0d0)**(i+j)*matout(i,j)
 20     continue
 10   continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccc
