cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     utils.f
c     RGF
c     modular version
c     Revision 2, Jul 5, 2001
c     Revision 2.2 for TASC:   Jul 23, 2001
c
c     Michael Hirsch / Wolfgang Quapp
c     Uni. Leipzig

      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

      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

      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

      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
      
      subroutine vec_iinit(iarray,ndim,iwert)
c Set all elements of ARRA1 to the integer value of WERT
      integer  iarray,iwert,ndim,i
      dimension iarray(ndim)
      if(ndim .le. 0) return
      do i=1,ndim
      iarray(i)=iwert
      end do
      return
      end
      
      subroutine vec_oinit(oarray,ndim,owert)
c Set all elements of ARRA1 to the logical value of WERT
      integer  ndim,i
      logical oarray,owert
      dimension oarray(ndim)
      if(ndim .le. 0) return
      do i=1,ndim
      oarray(i)=owert
      end do
      return
      end
      
      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

      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

      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
      
      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
      
      subroutine vec_mat_rentry(vec,dmat,npos,ndim,nrow,zkind)
c Copies VEC 'i'n or 'o'ut the NPOS row of matrix MAT     
      real*8 vec,dmat
      integer npos,ndim,nrow,i
      character zkind
      dimension vec(ndim),dmat(nrow,ndim)
      if(npos.gt.nrow) return
      if (zkind.eq.'i') then
      do 10,i=1,ndim
   10 dmat(npos,i)=vec(i)
      endif
      if (zkind.eq.'o') then
      do 20,i=1,ndim
   20 vec(i)=dmat(npos,i)
      endif
      return
      end

      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

      real*8 function adet(mat,sdet,dim)
c  calculate the Determinant of MAT      
      real*8 mat,amat,det,fact,sdet,sz
      integer dim,itest,sp,nexp,k,dreieck,ll
      dimension mat(dim,dim),amat(dim,dim),sp(dim-1),fact(dim)
      sz=1.0d0/16.0d0
      call mat_dcopy(mat,amat,dim,dim)
      itest=dreieck(amat,sp,dim)
      do 5,k=1,dim-1
   5  continue
      if (itest.eq.-1) then
      	adet=0.0d0
	sdet=0.0d0
      else
      	nexp=0
      	det=1
	sdet=1.0d0
      	do 10,k=1,dim
      		det=det*amat(k,k)
		fact(k)=amat(k,k)
      		if (k.lt.dim) then
      			if (k.lt.sp(k)) det=-det
      		endif
  20   		continue
      		if (dabs(det).lt.sz) then
      			det=det*16.0d0
      			nexp=nexp-1
      			goto 20
      		endif
  30   		continue
  		if (dabs(det).ge.1.0d0) then
  			det=det*sz
  			nexp=nexp+1
  			goto 30
      		endif
		if (k.eq.dim-1) then
			sdet=det
			ll=nexp
		endif
  10   	continue
      	if (nexp.lt.0) then
      		nexp=-nexp
      		do 40,k=1,nexp
      			det=det*sz
  40    	continue
  	elseif (nexp.gt.0) then 
  		do 50,k=1,nexp
  			det=det*16.0d0
  50		continue	
  	endif
	if (ll.lt.0) then
      		ll=-ll
      		do 410,k=1,ll
      			sdet=sdet*sz
  410    	continue
  	elseif (ll.gt.0) then 
  		do 510,k=1,ll
  			sdet=sdet*16.0d0
  510		continue	
  	endif
  	adet=det
      endif
      return
      end

      subroutine recorrect(corr,cold,gmat,ocorr,nn)
c   manipulation of the corrector for a better convergence
      integer nn
      real*8 corr,cold,dmnorm,gmat,vec_mdot
      logical olen,odir,ocorr,odone
      dimension corr(nn),cold(nn),gmat(nn,nn)
      ocorr=.false.
      if (dmnorm(corr,gmat,nn).gt.dmnorm(cold,gmat,nn)) then
           olen=.true.
      else 
           olen=.false.
      endif
      if (vec_mdot(corr,gmat,cold,nn).lt.0.0d0) then
           odir=.true.
      else 
           odir=.false.
      endif
      odone=.false.
      if (odir.and.olen.and.ocorr) then
        call vec_dmult_constant(corr,nn,0.6d0,corr)
	odone=.true.
      endif
      ocorr=odone
      return
      end
      

      integer function dreieck(mat,sp,dim)
c     Dreiecksfaktorisierung of MAT ; dim(SP) = DIM-1
      real*8 mat,z
      integer sp,dim,it,k,i,s,j
      dimension mat(dim,dim),sp(dim-1)
      it=0
      do 10,k=1,dim-1
      	z=0.0d0
      	do 20,i=k,dim
      		if (dabs(mat(i,k)).gt.z) then
      			z=dabs(mat(i,k))
      			s=i
      		endif
  20   	continue
  	if (z.eq.0.0d0) then
  		dreieck=-1
  		return
  	endif
  	sp(k)=s
  	if (k.lt.s) then
  		do 30,j=1,dim
  			z=mat(k,j)
  			mat(k,j)=mat(s,j)
  			mat(s,j)=z
  30		continue
  	endif
	do 40,i=k+1,dim
		mat(i,k)=mat(i,k)/mat(k,k)  
  40	continue
  	do 50,j=k+1,dim
  		do 60,i=k+1,dim
  			mat(i,j)=mat(i,j)-mat(i,k)*mat(k,j)
  60		continue
  50	continue
  10  continue    
      if (mat(dim,dim).eq.0.0d0) it=-1
      dreieck=it
      return
      end

      subroutine matadd(dm1,fact,dm2,out,n,m)
c    OUT(N,M) = DM1(N,M) + FACT * DM1(N,M)
      real*8 dm1,fact,dm2,out
      integer m,n
      dimension dm1(n,m),dm2(n,m),out(n,m)
      call vec_dmult_add(dm1,dm2,n*m,fact,out)
      return
      end

      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


       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          

      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
cc      integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
cc      common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
      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
      det=dmgs(basis,dmet,ndim)	
      call vec_dinit(compl,(ndim-1)*ndim,0.0d0)
      call mat_trans(basis,ndim,ndim,basis)
	do 160,i=2,ndim
      call vec_mrowcopy(i,basis,ndim,i-1,compl,ndim-1,ndim)
  160 continue
      return
      end
      
      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
      
      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/1.0d-55/
      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
      call cerr(-23)
      return
      end
      
      subroutine vec_madd
     c  (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

      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)
      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



      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)
      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
	       
	       
      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 50,m=k,n          
               sk=sk+w(m)*mataus(m,l)           
 50         continue          
            do 60,m=k,n          
               mataus(m,l)=mataus(m,l)-2*sk*w(m)          
 60         continue                         
 40       continue                          
 10   continue          
      return           
      end                                       

      subroutine closys
c closes internal files of mod. RGF: "store.rgf" "param.rgf" "messages.rgf"
      integer iochk
      integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
      common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
      close(isto,iostat=iochk)
      close(ipar,iostat=iochk)
      close(imsg,iostat=iochk)
      return
      end


      subroutine print_switch(io)
c   Output of the switches        
      integer io
      logical o_flat,      o_free, o_msg,    o_no_up,
     c        o_const_sgn, o_even, o_recorr, o_nn
      character*40 zkey
      common/switch/o_flat,o_free,o_msg,o_no_up,o_const_sgn,
     c  o_even,o_recorr,o_nn
      call keytable(zkey,14)
      if(o_flat)then
                write(io,10)zkey
      else
                write(io,20)zkey
      endif
      call keytable(zkey,15)
      if(o_free)then
                write(io,10)zkey
      else
                write(io,20)zkey
      endif
      call keytable(zkey,16)
      if(o_msg)then
                write(io,10)zkey
      else
                write(io,20)zkey
      endif
      call keytable(zkey,17)
      if(o_no_up)then
                write(io,10)zkey
      else
                write(io,20)zkey
      endif
      call keytable(zkey,18)
      if(o_const_sgn)then
                write(io,10)zkey
      else
                write(io,20)zkey
      endif
      call keytable(zkey,19)
      if(o_even)then
                write(io,10)zkey
      else
                write(io,20)zkey
      endif
      call keytable(zkey,20)
      if(o_recorr)then
                write(io,10)zkey
      else
                write(io,20)zkey
      endif
      call keytable(zkey,21)
      if(o_nn)then
                write(io,10)zkey
      else
                write(io,20)zkey
      endif
      return
  10  format(1x,A25,': On') 
  20  format(1x,A25,': Off') 
      end

      subroutine toang(vec1,vec2,nn)
c   (bohr,radian)->(angstroem,degree) for internal coordinates      
      real*8 vec1,vec2,torad,btoa
      integer i,nr,nn
      dimension vec1(nn),vec2(nn)
      data btoa  /0.52917706d0/
      torad=dacos(-1.0d0)/180.0d0
      nr=(nn+6)/3-1
      do i=1,nn
      if (i.le.nr) then
         vec2(i)=vec1(i)*btoa
      else 
         vec2(i)=vec1(i)/torad
      endif
      end do
      return
      end
      
      subroutine message(zstr,lstr)
c   writes the message ZSTR with length LSTR into imsg ("message.rgf")     
      character*40 zstr
      integer lstr
      integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
      logical o_flat,o_free,o_msg,o_no_up,o_const_sgn,o_even,o_recorr,
     c        o_nn
      common /switch/o_flat,o_free,o_msg,o_no_up,o_const_sgn,
     c               o_even,o_recorr,o_nn
      common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
      if (o_msg) write(imsg,zstr(1:lstr))
      return
      end
      
      subroutine err(iochk,zsource,lenz)
c  errorhandling       
      integer iochk,lenz
      integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
      character*40 zsource
      logical o_flat,o_free,o_msg,o_no_up,o_const_sgn,
     c               o_even,o_recorr,o_nn
      common /switch/o_flat,o_free,o_msg,o_no_up,o_const_sgn,
     c               o_even,o_recorr,o_nn
      common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
      if (iochk.ne.0) then
	 if (o_msg) then
             write(imsg,10)zsource(1:lenz)
	 else
             write(ires,10)zsource(1:lenz)
	 endif
      endif
      call cerr(iochk)
      return
  10  format(1x,'Error-Source: ',A)
      end
      
      subroutine cerr(iochk)
c  errorhandling       
      character*40 zkey
      integer iochk
      integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
      logical o_flat,o_free,o_msg,o_no_up,o_const_sgn,
     c               o_even,o_recorr,o_nn
      common /switch/o_flat,o_free,o_msg,o_no_up,o_const_sgn,
     c               o_even,o_recorr,o_nn
      common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
c     errorhandling: Number of error - IOCHK, outputfile - IOFILE 
      if ((.not.o_msg).and.(iochk.ne.0)) call exit (1)
      select case (iochk)
      case (0) 
         return
      case (1:40)
         write(imsg,10)iochk
      case (-1)
         write(imsg,20)
      case (-29:-10)
         call errtable(zkey,iochk)
         write(imsg,30)zkey
      case (-40:-30)
         call errtable(zkey,iochk)
         write(imsg,50)zkey
      case default
         write(imsg,40)iochk
      end select
      write(imsg,100)
      call closys
      call closext
      call exit(1)
      return       
  10  format(1x,'IO Error with number ',i3,' occured')
  20  format(1x,'Unexpected End-Of_File occoured')
  30  format(1x,'IO Error occured. Source: ',A40)
  40  format(1x,'strange IO Error occured. number:',i4)
  50  format(1x,'Internal numeric Error occured:',/,A40)
 100  format(1x,'*** Exit program with status 1***')
      end
       
      subroutine openext
c   check and open files for output, point, gradinet, hessian, b-matrix,      
c   and energy 
      integer iochk
      integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
      character*40 z_store,z_para,z_msg,z_FoR,z_FoP,z_FoG,z_FoH,z_FoB,
     c     z_iostat,z_FoE
      logical o_flat,opened,o_msg,o_free,o_no_up,o_const_sgn,
     c               o_even,o_recorr,o_nn
      common /switch/o_flat,o_free,o_msg,o_no_up,o_const_sgn,
     c               o_even,o_recorr,o_nn
      common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
      common/iofiles/z_store,z_para,z_msg,z_FoR,z_FoP,
     c              z_FoG,z_FoH,z_FoB,z_FoE
      opened=.false.
      if (o_msg) then
        inquire(file=z_msg,opened=opened)
        if (.not.opened) 
     c      open(imsg,file=z_msg,status='unknown',access='append')
        opened=.true.
      endif
      open (ires,file=z_FoR,status='unknown',access='append',
     c      err=10,iostat=iochk)
      if (opened) write(imsg,*)z_FoR//z_iostat(iochk)
      open (ipnt,file=z_FoP,err=10,status='old',iostat=iochk)
      if (opened) write(imsg,*)z_FoP//z_iostat(iochk)
      open (igrd,file=z_FoG,err=10,status='old',iostat=iochk)
      if (opened) write(imsg,*)z_FoG//z_iostat(iochk)
      open (ihes,file=z_FoH,err=10,status='old',iostat=iochk)
      if (opened) write(imsg,*)z_FoH//z_iostat(iochk)
      if (.not.o_flat) then
        open (iwil,file=z_FoB,err=10,status='old',iostat=iochk)
        if (opened) write(imsg,*)z_FoB//z_iostat(iochk)
      endif
      open (iene,file=z_FoE,err=10,status='old',iostat=iochk)
      if (opened) write(imsg,*)z_FoE//z_iostat(iochk)

      return
   10 continue
      call cerr(-12)
      return
      end		
      
      subroutine closext
c   close files for output, point, gradinet, hessian, b-matrix      
      integer iochk
      integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
      common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene
      close(ires,iostat=iochk)
      close(ipnt,iostat=iochk)
      close(igrd,iostat=iochk)
      close(ihes,iostat=iochk)
      close(iwil,iostat=iochk)
      close(iene,iostat=iochk)
      return
      end
        
      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
      
      real*8 function diag_inv(diag,dinv,nn)
c  inverse DINV of a diagonalmatrix DIAG      
      real*8 diag,dinv,det
      integer nn,i
      dimension diag(nn,nn),dinv(nn,nn)
      call vec_dinit(dinv,nn*nn,0.0d0)
      det=1.0d0
      do i=1,nn
      if (diag(i,i).eq.(0.0d0)) call cerr(-30)
      det=det*diag(i,i)
      dinv(i,i)=1.0d0/diag(i,i)
      enddo
      diag_inv=det
      return
      end
      
      logical function omat_inv(in,ndim,out)
c    calculates the inverse matrix of IN       
      real*8 in,out
      integer ndim,i,j,k
      dimension in(ndim,ndim),out(ndim,ndim)
      call mat_dcopy(in,out,ndim,ndim)
      do 10,k=1,ndim
         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,ndim
	     if (i.ne.k) out(i,k)=-(out(i,k)*out(k,k))
   20    continue
         do 30,j=1,ndim
	     if (j.ne.k) then
	        do 40,i=1,ndim
		      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

      subroutine getmet(bmat,bpinv,gmat,ginv,nn)
c  calculation of GMAT and GINV      
      real*8 bmat,bpinv,gmat,ginv
      integer nn
      logical omat_inv
      dimension bmat(nn,nn+6),gmat(nn,nn),ginv(nn,nn),bpinv(nn,nn+6)
c  build metric, its inverse, the pseudo inverse of B-matrix from B-matrix
      call vec_dinit(gmat,nn*nn,0.0d0)
      call blas_dgemm('n','t',nn,nn,nn+6,1.0d0,bmat,nn,
     c     bmat,nn,1.0d0,gmat,nn)
      if (.not.omat_inv(gmat,nn,ginv)) call err(-32,'getmet',6)  
      call vec_dinit(bpinv,nn*(nn+6),0.0d0)  
      call blas_dgemm('n','n',nn,nn+6,nn,1.0d0  
     c     ,ginv,nn,bmat,nn,1.0d0,bpinv,nn)   
      return
      end
      
      subroutine blas_dgemm
c    matrix multiplikation      
     *     (mata,matb,nrow,ncol,mdi,zahl1,arra1,
     *     ndima,arra2,ndimb,zahl2,arra3,ndimc)
      character*1 mata,matb
      integer*4   ncol,nrow,ndima,ndimb,ndimc,mdi
      real*8      zahl1,zahl2,arra1
      real*8      arra2,arra3
      real*8      zwischen
      integer*4  i,j,k
      dimension arra1(*),arra2(*),arra3(*)  
      if (mdi.eq.0) then
         do 1,i=1,nrow
            do 2,j=1,ncol
               arra3(i+(j-1)*ndimc)=zahl2*arra3(i+(j-1)*ndimc)
 2          continue
 1       continue
      return
      endif
      if (mata.eq.'n'.and.matb.eq.'n') then
          do 3,i=1,nrow
             do 4,j=1,ncol
                zwischen = 0.0
                do 5,k=1,mdi
                   zwischen =  arra1(i+(k-1)*ndima) *
     c             arra2(k+(j-1)*ndimb) + zwischen
 5              continue
                arra3(i+(j-1)*ndimc) = zahl1*zwischen + 
     c           zahl2*arra3(i+(j-1)*ndimc)
 4           continue
 3        continue
          return
      endif
      if (mata.eq.'n'.and.matb.eq.'t') then
          do 6,i=1,nrow
             do 7,j=1,ncol
                zwischen = 0.0
                do 8,k=1,mdi
                   zwischen =  arra1(i+(k-1)*ndima) *
     c             arra2(j+(k-1)*ndimb) + zwischen
 8              continue
                arra3(i+(j-1)*ndimc) = zahl1*zwischen +
     c           zahl2*arra3(i+(j-1)*ndimc)
 7           continue
 6        continue
          return
      endif
      if (mata.eq.'t'.and.matb.eq.'n') then
          do 9,i=1,nrow
             do 10,j=1,ncol
                zwischen = 0.0
                do 11,k=1,mdi
                   zwischen =  arra1(k+(i-1)*ndima) *
     c             arra2(k+(j-1)*ndimb) + zwischen
 11              continue
                arra3(i+(j-1)*ndimc) = zahl1*zwischen +
     c           zahl2*arra3(i+(j-1)*ndimc)
 10          continue
 9        continue
          return
      endif
      if (mata.eq.'t'.and.matb.eq.'t') then
          do 12,i=1,nrow
             do 13,j=1,ncol
                zwischen = 0.0
                do 14,k=1,mdi
                   zwischen =  arra1(k+(i-1)*ndima) *
     c             arra2(j+(k-1)*ndimb) + zwischen
 14              continue
                arra3(i+(j-1)*ndimc) = zahl1*zwischen +
     c           zahl2*arra3(i+(j-1)*ndimc)
 13          continue
 12       continue
          return
      endif
      print*,'blas_dgemm failed'
      stop
      return
      end


      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
      
      subroutine nka(ginv,hi,eival,eivec,nn)
c   normal coordinate analysis      
      real*8 schol,ginv,hi,tmp1,tmp2,eival,eivec
      integer nn
      dimension ginv(nn,nn),schol(nn,nn),tmp1(nn,nn),tmp2(nn,nn)
      dimension eival(nn),eivec(nn,nn)
      call vec_dcopy(ginv,schol,nn*nn)  
      call cholesky(schol,nn)  
      call vec_dinit(tmp1,nn*nn,0.0d0)  
      call blas_dgemm  
     $     ('n','n',nn,nn,nn,1.0d0,hi,  
     $     nn,schol,nn,1.0d0,tmp1,nn)  
      call vec_dinit(tmp2,nn*nn,0.0d0)  
      call blas_dgemm  
     $     ('t','n',nn,nn,nn,1.0d0,schol,  
     $     nn,tmp1,nn,1.0d0,tmp2,nn)  
      call eigen(tmp2,nn,nn,eival,eivec,tmp2,0)  
      call matmult(schol,nn,eivec,nn,eivec,nn)
      return
      end
      
      integer function i_dfp(einmat,q,p,nn,upmat)
c   DFP-Update of Hessian, see     
c Heidrich,Kliesch,Quapp: Properties of Potential Energy Surfaces...
c S.51 Eqn.(21)/(22)
c
c einmat	Hessian to update      (H_0)
c q		difference of gradient (g_1 - g_0)
c p		difference of point    (x_1 - x_0)
c upmat         updated Hessian        (H_1)
      real*8 einmat,q,p,upmat,qp,unit,b1,b2,b3,
     1     vec_ddot,eps,sw
      integer nn
      dimension unit(nn,nn),b3(nn,nn),
     +  einmat(nn,nn),upmat(nn,nn),
     +  b1(nn,nn),b2(nn,nn),sw(nn,nn),
     +  q(nn), p(nn)
      data eps/1.0d-30/
      call vec_dinit(upmat,nn*nn,0.0d0)
      call mat_diag(unit,nn,1.0d0)
      qp=vec_ddot(q,p,nn)
      if(dabs(qp).lt.eps) call cerr(-31)
      call matmult(q,nn,p,nn,sw,1)
      call matadd(unit,(-1.0d0)/qp,sw,b1,nn,nn)
      call matmult(p,nn,q,nn,sw,1)
      call matadd(unit,(-1.0d0)/qp,sw,b2,nn,nn)
      call matmult(q,nn,q,nn,b3,1)
      call tri_mat(b1,nn,nn,einmat,nn,b2,nn,sw)
      call matadd(sw,1.0d0/qp,b3,upmat,nn,nn)
      i_dfp=0
    2 continue
      return
      end
      
      real*8 function newton(hi,fi,step,gmat,nn)
c  calculates STEP = - HI^(-1) * FI and returns the steplength
      real*8 hi,fi,step,gmat,gl,vec_mdot
      integer i,nn
      dimension hi(nn,nn),fi(nn),step(nn),gmat(nn,nn),gl(nn,nn+1)
      call mat_dcopy(hi,gl,nn,nn)
      do 10,i=1,nn
         gl(i,nn+1)=-fi(i)
 10   continue
      call linsolve(gl,step,nn)       
      newton=dsqrt(vec_mdot(step,gmat,step,nn))     
      return
      end
      
      subroutine sig(ev,iv,nn)
c count the number IV of positive, negative and zero eigenvalues EV
      real*8 ev
      integer iv,nn,i
      dimension ev(nn),iv(3)
      do i=1,3
      iv(i)=0
      enddo
      do i=1,nn
      if (ev(i).gt.(0.0d0)) iv(1)=iv(1)+1
      if (ev(i).lt.(0.0d0)) iv(2)=iv(2)+1
      if (ev(i).eq.(0.0d0)) iv(3)=iv(3)+1
      enddo
      return
      end
      
      subroutine keytable(zkey,nr)
c   Table of KEYNAMES for the parameterfile      
      character*40 zkey,ztab
      integer nr
      dimension ztab(30)
      ztab(1) ='Direction_of_Search'
      ztab(2) ='Dimension'
      ztab(3) ='Job_Name'
      ztab(4) ='File_of_Points'
      ztab(5) ='File_of_Gradient'
      ztab(6) ='File_of_Hessian'
      ztab(7) ='File_of_Results'
      ztab(8) ='File_of_B_matrix'
      ztab(9) ='Stopping_Criterion'
      ztab(10)='Predictor_Corrector_Criterion'
      ztab(11)='Length_of_Predictorsteps'
      ztab(12)='Min_Number_of_Steps'
      ztab(13)='Max_Number_of_Steps'
      ztab(14)='Flat_Space'
      ztab(15)='Free_First_Point'
      ztab(16)='Messages'
      ztab(17)='No_Update'
      ztab(18)='Const_Sign_Update'
      ztab(19)='Search_Even_Index'
      ztab(20)='Recorrect'
      ztab(21)='nul_switch'
      ztab(22)='nul'
      ztab(23)='nul'
      ztab(24)='nul'
      ztab(25)='nul'
      ztab(26)='nul'
      ztab(27)='nul'
      ztab(28)='nul'
      ztab(29)='File_of_Energy'
      ztab(30)='nul'
      zkey=ztab(nr)
      return
      end
      
      
      subroutine errtable(zkey,nr)
c     table for errormessages      
      character*40 zkey,ztab
      integer nr,i
      dimension ztab(40)
   10 continue  
c               1        0         0         0          0                
      ztab(1)='unspecific'
      ztab(2)='unspecific'
      ztab(3)='unspecific'
      ztab(4)='unspecific'
      ztab(5)='unspecific'
      ztab(6)='unspecific'
      ztab(7)='unspecific'
      ztab(8)='unspecific'
      ztab(9)='unspecific'
      ztab(10)='Direction_of_Search'
      ztab(11)='Number_of_Atoms'
      ztab(12)='Job_Name'
      ztab(13)='File_of_Points'
      ztab(14)='File_of_Gradient'
      ztab(15)='File_of_Hessian'
      ztab(16)='File_of_Results'
      ztab(17)='File_of_B_matrix'
      ztab(18)='Stopping_Criterion'
      ztab(19)='Predictor_Corrector_Criterion'
      ztab(20)='Length_of_Predictorsteps'
      ztab(21)='Min_Number_of_Steps'
      ztab(22)='Max_Number_of_Steps'
      ztab(23)='Gram-Schmidt algorithm'
      ztab(24)='unspecific'
      ztab(25)='unspecific'
      ztab(26)='unspecific'
      ztab(27)='unspecific'
      ztab(28)='unspecific'
      ztab(29)='File_of_Energy'
      ztab(30)='Div. by zero while inverting diag-matrix'
      ztab(31)='DFP update failed'
      ztab(32)='matrix inverting fails'
      ztab(33)='file not exists'
      ztab(34)='unspecific'
      ztab(35)='unspecific'
      ztab(36)='unspecific'
      ztab(37)='unspecific'
      ztab(38)='unspecific'
      ztab(39)='unspecific'
      ztab(40)='unspecific'
      if (nr.eq.0) then
        i=10
	do while ((zkey.ne.ztab(i)).and.(i.le.22))
	i=i+1
	end do
	nr=-i
	return
      endif
      if ((nr.lt.0).and.(nr.ge.-40)) then
        zkey=ztab(-nr)
	return
      endif
      zkey='not known'
      return
      end

      
      
       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          

       subroutine cholesky(g,ndim)  
c  cholesky decomposition       
       real*8 g,d  
       integer ndim,i,j,k,l  
       dimension g(ndim,ndim)  
       i=ndim-1  
       do 1,j=1,i  
          k=j+1  
          do 2,l=k,ndim  
             g(j,l)=0.0d0  
 2        continue  
 1     continue  
       do 3,i=1,ndim  
          if (g(i,i).le.0.0d0) print*,'ERROR CHOL'  
          g(i,i)=dsqrt(g(i,i))  
          j=i+1  
          if (j.gt.ndim) goto 3  
          d=g(i,i)  
          do 4,k=j,ndim  
             g(k,i)=g(k,i)/d  
 4        continue  
          do 5,k=j,ndim  
             do 6,l=k,ndim  
                g(l,k)=g(l,k)-g(l,i)*g(k,i)  
 6           continue  
 5        continue  
 3     continue       
       return  
       end   
      
      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


      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)          
      call loesung(dre,x,n)          
      return          
      end          
  
       subroutine qr(gex,q,n)        
C      QR-Zerlegung der Jacobimatrix von ge (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 
       return        
       end  
  
  
      subroutine strcp(z1,z2,nn)
c     copies Character*NN   Z1 to Z2      
      character*1 z1,z2
      integer nn,i
      dimension z1(nn),z2(nn)
      do 10,i=1,nn
   10 z2(i)=z1(i)
      return
      end
      
      subroutine recreate(ifile,zfile)
c     create an empty file
      integer iochk,ifile
      character*40 zfile
      open(ifile,file=zfile,iostat=iochk,status='unknown')
      close(ifile,iostat=iochk,status='delete')
      return
      end      

      subroutine getvec(iofile,nkey,vec,nn)
c reads vector VEC following the key with number NKEY from IOFILE     
      logical oready
      integer nkey,i,iochk,nn,iofile
      character*40 zkey,zscr
      real*8 vec
      dimension vec(nn)
      oready=.false.
      rewind(iofile)
      call keytable(zscr,nkey)
      do while (.not.oready)
      read(iofile,*,end=10,iostat=iochk)zkey
      if (zkey.eq.zscr) then
          i=1
	  do while (i.le.nn)
	    read(iofile,*,iostat=iochk)vec(i)
	    call cerr(iochk)
	    i=i+1
	  end do
	  oready=.true.
      endif
      end do
  10  continue
      if (.not.oready) call cerr(-10)
      return
      end
      
      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
      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
     c   (tangent,nn,1/dmnorm(tangent,gmat,nn),tangent)
      call mat_dcopy(curveq,ss1,nn-1,nn)
      do 10,i=1,nn            
   10    ss1(nn,i)=tangent(i)         
      if (det(ss1,nn).lt.zero) then
        call vec_dmult_constant(tangent,nn,-1.0d0,tangent)
      endif
      return
      end      
      
      subroutine printv(io,val,zval,m)
c   writes a table ZVAL, VAL with M rows in IO 
      integer m,io,i
      real*8 val
      character*20 zval
      dimension val(m),zval(m)
      write(io,20)
      write(io,30)
      write(io,20)
      do i=1,m
      write(io,10)zval(i),val(i)
      end do
      write(io,20)
      return
 10   format (1x,A20,1x,f8.5)
 20   format (1x,29('-'))
 30   format (12x,'scalars')
      end
      
      real*8 function snorm(adj,grad,ginv,gmat,nn)
c   S norm      
      real*8 adj,grad,ginv,dmnorm,ag,gmat
      integer nn
      dimension adj(nn,nn),grad(nn),ginv(nn,nn),ag(nn),gmat(nn,nn)
      call matmult(adj,nn,grad,1,ag,nn)
      snorm=dmnorm(ag,gmat,nn)/dmnorm(grad,ginv,nn)
      return
      end
      
      
      real*8 function bnorm(adj,grad,gmat,ag,nn)
c   Branin norm      
      real*8 adj,grad,gmat,ag,dmnorm
      integer nn
      dimension adj(nn,nn),grad(nn),gmat(nn,nn),ag(nn)
      call matmult(adj,nn,grad,1,ag,nn)
      bnorm=dmnorm(ag,gmat,nn)
      call vec_dmult_constant(ag,nn,-1.0d0,ag)
      return
      end

      real*8 function ray(adj,grad,ginv,nn)
c rayleigh coef.      
      real*8 adj,grad,ginv,ag,vec_ddot,vec_mdot
      integer nn
      dimension adj(nn,nn),grad(nn),ginv(nn,nn),ag(nn)
      call matmult(adj,nn,grad,1,ag,nn)
      ray=vec_ddot(grad,ag,nn)/vec_mdot(grad,ginv,grad,nn)
      return
      end
      
      subroutine mat_adj(matin,ndim,matout)
c adjoint matrix      
      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

      
      real*8 function cnorm(vec,nn)
c   cartesian NORM of VEC      
      real*8 vec,vec_ddot
      integer nn
      dimension vec(nn)
      cnorm=dsqrt(vec_ddot(vec,vec,nn))
      return
      end  

      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  
      
      subroutine help(io,index)
      integer io,index
c     5=stdin, 6=stdout, 0=stderror      
      character*120 zpara,zhelp,zver,zret
      character*20 zrgf
      dimension zpara(5),zhelp(3),zver(5),zret(5)

      zpara(1)=
     c'Format of param.rgf:\n (Components of searchdir separated '//
     c'by space or newline)\n &param\n job_id=(number)'
      zpara(2)=
     c'nn=(dimension)\n z_FoP=(filename)\n z_FoG=(filename)\n z_FoH'//
     c'=(filename)\n z_FoB=(filename)\n z_res=(filename)'
      zpara(3)=
     c 'o_flat={T|F}\n o_free={T|F}\n o_msg={T|F}\n o_no_up={T|F}\n'//
     c ' o_const_sgn={T|F}\n o_even={T|F}\n o_recorr={T|F}'
      zpara(4)=
     c 'o_nn={T|F}\n crit_stop=(Fortran double)\n crit_pc=(Fortran'//
     c ' double)\n pstepl=(Fortran double)'
      zpara(5)=
     c 'itmin=(stepnumber)\n itmax=(stepnumber)\n /' 

      zhelp(1)='Usage: (programname) [options]\n Options:\n '//
     c   '-ver\t\tprint version information\n '//
     c   '-para\t\tprint the format of »param.rgf«'
      zhelp(2)='-ret\t\tprint a desciption of return values\n '//
     c '-exact\t\tforces "exact HESSIAN"\n '//
     c '-update\tforces "update HESSIAN"'
      zhelp(3)='-reinit\treread PSTEPL,CRIT_STOP,CRIT_PC,ITMIN,ITMAX'//
     & '\n\t\tfrom the parameterfile'
     

      zver(1)='RGF (»Reduced Gradient Following«) is a realisation'//
     c '\n of the "following a reduced gradient" - method\n'
      zver(2)=
     c 'Reference:\tW.Quapp, M.Hirsch, D.Heidrich:\n '//
     c ' \t\tTheor Chem Acc (1998) 100:285-299'
      zver(3)=
     c '\t\tW.Quapp, M.Hirsch, O.Imig, D.Heidrich:\n \t\tJ Comput '//
     c 'Chem 19 (1998) 1087-1100\n'
      zver(4)=
     c 'First release: Apr 19, 2000\n Revision 1:    Nov 11, 2000\n'//
     c ' Revision 2:    Jul 05, 2001\n'
      zver(5)='Contact to the authors and programmers:\n '//
     c 'E-mail:\thirsch@mathematik.uni-leipzig.de\n '//
     c '\t\tquapp@mathematik.uni-leipzig.de' 
     

      zret(1)=
     c 'Legend of return values\n 0\tafter printing any help message'//
     c ', like this\n 1\tparam.rgf not exists'
      zret(2)=
     c '\tOR error while reading parameters\n 4\tPREDICTOR in step '//
     c 'executed;\n\tnext hessian: UPDATE'
      zret(3)=
     c '5\tPREDICTOR in step executed;\n\tnext hessian: EXACT\n '//
     c '6\tCORRECTOR in step executed;\n\tnext hessian: UPDATE'
      zret(4)=
     c '7\tCORRECTOR in step executed;\n\tnext hessian: EXACT \n'//
     c '10\tmax. Stepnumber reached'
      zret(5)=
     c '\tOR stopping criterion satisfied'


      call getarg (0,zrgf)

      if (index.eq.0) write(io,*) zhelp
      if (index.eq.1) write(io,*) zpara
      if (index.eq.2) write(io,*) zver
      if (index.eq.3) write(io,*) zret
      return
      end
                        
      character*40 function z_iostat(iostat)		
      integer iostat
      select case (iostat)
      case(0)  
                 z_iostat='successfull io'
      case(100)
                 z_iostat='error in format'
      case(101)
       z_iostat='illegal unit number'
      case(102)
       z_iostat='formatted io not allowed'
      case(103)
       z_iostat='unformatted io not allowed'
      case(104)
       z_iostat='direct io not allowed'
      case(105)
       z_iostat='sequential io not allowed'
      case(106)
       z_iostat='cant backspace file'
      case(107)
       z_iostat='null file name'
      case(108)
       z_iostat='cant stat file'
      case(109)
       z_iostat='unit not connected'
      case(110)
       z_iostat='off end of record'
      case(111)
       z_iostat='truncation failed in endfile'
      case(112)
       z_iostat='incomprehensible list input'
      case(113)
       z_iostat='out of free space'
      case(114)
       z_iostat='unit not connected'
      case(115)
       z_iostat='read unexpected character'
      case(116)
       z_iostat='bad logical input field'
      case(117)
       z_iostat='bad variable type'
      case(118)
       z_iostat='bad namelist name'
      case(119)
       z_iostat='variable not in namelist'
      case(120)
       z_iostat='no end record'
      case(121)
       z_iostat='variable count incorrect'
      case(122)
       z_iostat='subscript for scalar variable'
      case(123)
       z_iostat='invalid array section'
      case(124)
       z_iostat='substring out of bounds'
      case(125)
       z_iostat='subscript out of bounds'
      case(126)
       z_iostat='cant read file'
      case(127)
       z_iostat='cant write file'
      case(128)
       z_iostat='\"new\" file exists'
      case(129)
       z_iostat='cant append to file'
      case(130)
       z_iostat='non-positive record number'
      case(131)
       z_iostat='I/O started while already doing I/O'
      case default
       z_iostat='unknown IO status'
      end select
      return
      end

      subroutine handleargs(n,ostat)
      integer n,i,istat,ind
      logical ostat
      character*10 zarg
      character*40 zall
      dimension ostat(8)
      data zall /'-ver-para-ret-exact-update-reinit'/
      do i=3,8
         ostat(i)=.false.
      enddo	 
      do i=1,n
        call getarg(i,zarg)
	istat=len_trim(zarg)
	ind=index(zall,zarg(1:istat))
        select case (ind)
	  case(1)
	     call help(6,2)
	     call exit(0)
	  case(5)
	     call help(6,1)
	     call exit(0)
	  case(10)
	     call help(6,3)
	     call exit(0)
	  case(14)
	     ostat(3)=.true.
	  case(20)
	     ostat(4)=.true.
	  case(27)
	     ostat(5)=.true.
	  case default
	     call help(6,0)
	     call exit(0)
	end select
      enddo  
      return
      end

c     Modification since  
c     Thu Jun 19 13:01:45 CEST 2001

      subroutine func_prop(p1,p2,p3,iprop,zprop) 
      real*8 p1,p2,p3
      integer iprop,ibset
      character*20 zprop
      character*6 zmon,zext,zchg
      iprop=0
      zmon='decr  '
      zext='mono  '
      zchg='const '
      if (p3.gt.p2) then
          iprop=ibset(iprop,0)
	  zmon='incr  '
      endif	  
      if ((p2-p1)*(p2-p3).ge.0.0d0) then
          iprop=ibset(iprop,1)
	  zext='extr  '
      endif  
      if (p2*p3.lt.0.0d0) then
          iprop=ibset(iprop,2)
	  zchg='change'
      endif
      zprop=zmon//','//zext//','//zchg	  
      return
      end
      
      subroutine reread(zfile)
      integer job_id,nn,itmin,itmax,is
      character*40 z_FoP,z_FoG,z_FoH,z_FoB,z_FoR,zfile,z_FoE
      real*8 crit_stop,crit_pc,pstepl
      logical o_flat,o_free,o_msg,o_no_up,o_const_sgn,o_even,
     c o_recorr,o_nn,ostat
      dimension ostat(8)
      common /para/crit_stop,crit_pc,pstepl,itmin,itmax,ostat
      namelist /param/job_id,nn,z_FoP,z_FoG,z_FoH,z_FoB,z_FoR,
     c o_flat,o_free,o_msg,o_no_up,o_const_sgn,o_even,o_recorr,o_nn,
     c crit_stop,crit_pc,pstepl,itmin,itmax,z_FoE
      inquire(file=zfile,exist=ostat(1),opened=ostat(2),number=is)
      if (.not.ostat(1)) return
      if (ostat(2)) close(is)
      open(33,file=zfile,status='old')      
      read(33,param,iostat=is)
      if (is.ne.0.and.is.ne.119) then
         write(*,*)'wrong format of parameter file'
	 call help(6,1)
         call exit(1)
      endif
      close(33)
      return 
      end

C  SLSLSLSLSLSL SLSLSLSLSLSL SLSLSLSLSLSL SLSLSLSLSLSL SLSLSLSLSLSL
c
c  see http://csep1.phy.ornl.gov/CSEP/MO/NODE11B.html   or 
c      T Schlick, in Rev.Comput.Chem.,Vol.III (1992) p.1-72 
c
c SLSLSLSLSLSL SLSLSLSLSLSL SLSLSLSLSLSL SLSLSLSLSLSL SLSLSLSLSLSL
c
      subroutine search_line(o_even,gra,yy,en,fL2,po,xx,lam1,Stepl,
     c                       SeDir,nn)
       implicit real*8 (a-h,p-z), integer(i-n)
       real*8 lam1,gra,en,po,Stepl,SeDir,xx,yy,aL,bL,fL,fL2,gL,gL2,
     c        salt,vec_ddot,lambda
       integer nn 
       logical o_even
       dimension gra(nn),po(nn),SeDir(nn),xx(nn),yy(nn)
c
c      lam1 first EigenValue
c      xx line search point
c      yy gradient there
c      fL2 energy in test point
       salt=Stepl
       lambda=lam1
       fL= en
       gL= vec_ddot(gra,SeDir,nn)
       call vec_dmult_add(po,SeDir,nn,Stepl,xx) 
       gL2= vec_ddot(yy,SeDir,nn)
2      continue
      if(o_even) then 
       If( lambda .gt. 0.0d0 ) then
         aL= (-2.0d0*(fL2-fL)+(gL2+      gL)*Stepl)/Stepl**3
         bL= ( 3.0d0*(fL2-fL)-(gL2+2.0d0*gL)*Stepl)/Stepl**2
         If((aL.ne. 0.0d0).and.(bL**2-3.0d0*aL*gL.gt.0.0d0) ) then
          Stepl= -(bL +Dsqrt(bL**2-3.0d0*aL*gL))/(3.0d0*al)
         else
          Stepl= -gL/(2.0d0*bl)
         endif
       else
         If(fL2.gt. fL) then
         lambda=1.0d0
         goto 2
        else
         If(fL2.lt. fL) Stepl=Stepl*1.111111d0
        endif
       endif
c
      else   
       If(lambda.lt. 0.0d0) then
        aL=(-2.0d0*(fL2-fL)+(gL2+      gL)*Stepl)/Stepl**3
        bL=( 3.0d0*(fL2-fL)-(gL2+2.0d0*gL)*Stepl)/Stepl**2
         If((aL.ne. 0.0d0) .and. (bL**2-3.0d0*aL*gL.gt.0.0d0)) then
         Stepl= -(bL +Dsqrt(bL**2-3.0d0*aL*gL))/(3.0d0*al)
         else
         Stepl= -gL/(2.0d0*bl)
         endif
        else
cccccc
         If(fL2.gt.fL) Stepl=Stepl
c uphill in minimum bowl: LineSearch is still open !!
cccccc
        endif
      endif
c
       if(Stepl.gt.2.222d0*salt) Stepl=2.222d0*salt
       if(Stepl.lt.0.0d0)        Stepl=salt/2.0d0
       return
       end
