c    @(#)invmtx.f	1.2
c
      subroutine invmtx (a,na,v,nv,n,d,ip,ier)
      integer na,nv,n,ip(2*n),ier
      real a(na,n),v(nv,n),d
      data iexmax/75/
  115 format(28h0*matrix singular in invmtx*)
  116 format(34h0*determinant too large in invmtx*)
      ier = ierinv(n,na,nv)
      if (ier .ne. 0) return
      do 102 j=1,n
         ip(j) = 0
         do 101 i=1,n
            v(i,j) = a(i,j)
  101    continue
  102 continue
      d = 1.
      iex = 0
      do 110 m=1,n
         vmax = 0.
         do 104 j=1,n
            if (ip(j) .ne. 0) go to 104
            do 103 i=1,n
               if (ip(i) .ne. 0) go to 103
               vh = abs(v(i,j))
               if (vmax .gt. vh) go to 103
               vmax = vh
               k = i
               l = j
  103       continue
  104    continue
         ip(l) = k
         npm = n+m
         ip(npm) = l
         d = d*v(k,l)
  105    if (abs(d) .le. 1.0) go to 106
         d = d*0.1
         iex = iex+1
         go to 105
  106    continue
         pvt = v(k,l)
         if (m .eq. 1) pvtmx = abs(pvt)
         if (abs(pvt/float(m))+pvtmx .eq. pvtmx) go to 113
         v(k,l) = 1.
         do 107 j=1,n
            hold = v(k,j)
            v(k,j) = v(l,j)
            v(l,j) = hold/pvt
  107    continue
         do 109 i=1,n
            if (i .eq. l) go to 109
            hold = v(i,l)
            v(i,l) = 0.
            do 108 j=1,n
               v(i,j) = v(i,j)-v(l,j)*hold
  108       continue
  109    continue
  110 continue
      m = n+n+1
      do 112 j=1,n
         m = m-1
         l = ip(m)
         k = ip(l)
         if (k .eq. l) go to 112
         d = -d
         do 111 i=1,n
            hold = v(i,l)
            v(i,l) = v(i,k)
            v(i,k) = hold
  111    continue
  112 continue
      if (iex .gt. iexmax) go to 114
      d = d*10.**iex
      return
  113 ier = 33
      print 115
      return
  114 ier = 1
      d = float(iex)
      print 116
      return
      end
