      subroutine ginverse(n,gmat,ginvmat,wkmat,coeff,veci,vecj)

c  this subroutine written by alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      dimension gmat(n,*),ginvmat(n,*),wkmat(n,*),coeff(n,*)
      dimension veci(*),vecj(*)

      data loopmax/500/

      data big,tol,eigen/10000.0d0,0.00000001d0,0.001d0/
      data zero,half,one,two,hundred/0.0d0,0.5d0,1.0d0,2.0d0,100.0d0/

      do 1001 j=1,n
      do 1002 i=1,n
 1002 coeff(i,j)=zero
 1001 coeff(j,j)=one

      do 1003 j=1,n
      do 1003 i=1,n
      wkmat(i,j)=zero
      do 1003 k=1,n
 1003 wkmat(i,j)=wkmat(i,j)+coeff(k,j)*gmat(k,i)

      do 1004 loop=1,loopmax

      g_ij_max=zero

      do 1005 i=1,n
      do 1005 j=1,i

      if(i.eq.j) goto 1005

      g_ij=zero

      do 1006 k=1,n
 1006 g_ij=g_ij+coeff(k,j)*wkmat(k,i)

      if(abs(g_ij).gt.g_ij_max) g_ij_max=abs(g_ij)

      t01=tol/hundred/hundred/hundred

      if(abs(g_ij).gt.t01) then
                                 g_ii=zero

                                 do 1007 k=1,n
 1007                            g_ii=g_ii+coeff(k,i)*wkmat(k,i)

                                 g_jj=zero

                                 do 1008 k=1,n
 1008                            g_jj=g_jj+coeff(k,j)*wkmat(k,j)

                                 theta=(g_jj-g_ii)/(two*g_ij)
                           else
                                 goto 1005
                           endif

      if(theta.lt.zero) then
                              tan_phi=theta+sqrt(theta*theta+one)
                        else
                              tan_phi=theta-sqrt(theta*theta+one)
                        endif

      if(abs(theta).gt.big) tan_phi=-(half/theta)

      phi=atan(tan_phi)

      c=cos(phi)
      s=sin(phi)

      tau=s/(one+c)

      do 1009 k=1,n
 1009 veci(k)=coeff(k,i)

      do 1010 k=1,n
 1010 vecj(k)=coeff(k,j)

      do 1011 k=1,n
 1011 coeff(k,i)=coeff(k,i)+s*(vecj(k)-tau*veci(k))

      do 1012 k=1,n
 1012 coeff(k,j)=coeff(k,j)-s*(veci(k)+tau*vecj(k))

      do 1013 k=1,n
 1013 veci(k)=wkmat(k,i)

      do 1014 k=1,n
 1014 vecj(k)=wkmat(k,j)

      do 1015 k=1,n
 1015 wkmat(k,i)=wkmat(k,i)+s*(vecj(k)-tau*veci(k))

      do 1016 k=1,n
 1016 wkmat(k,j)=wkmat(k,j)-s*(veci(k)+tau*vecj(k))

 1005 continue

      if(g_ij_max.lt.tol) goto 1017

      if(loop.eq.loopmax) write(6,1018)
 1018 format(/,' ginverse routine did not fully converge')

      if(loop.eq.loopmax) stop

 1004 continue

 1017 continue

      do 1019 j=1,n
      do 1019 i=1,n
 1019 ginvmat(i,j)=zero

      iswitch=0

      do 1020 k=1,n

      g_value=zero

      do 1021 j=1,n
 1021 g_value=g_value+coeff(j,k)*wkmat(j,k)

      if(g_value.lt.eigen) goto 1020

      factor=sqrt(one/g_value)

      do 1022 j=1,n
 1022 coeff(j,k)=coeff(j,k)*factor

      do 1023 j=1,n
      do 1023 i=1,n
 1023 ginvmat(i,j)=ginvmat(i,j)+coeff(i,k)*coeff(j,k)

 1020 continue

      return
      end
