       subroutine invmtx2(a,n,nobs)
c ***
c *** subroutine invmtx2(a,n,nobs)
c *** invert the (nobs x nobs) matrix a, which is store in a nxn array.
c *** Use LU decomposition and an interative improvement
c ***/
c example driver:
c      program xinvmtx2
       include 'parameter.const'
c
c       parameter (n=5,nobs=5)
c      integer indx(nobs)
c      real d a(n,n)
c      DATA a/1.0,2.0,1.0,4.0,5.0,2.0,3.0,1.0,5.0,1.0,
c     *     3.0,4.0,1.0,1.0,2.0,4.0,5.0,1.0,2.0,3.0,
c     *     5.0,1.0,1.0,3.0,4.0/
c 101  format(1x,5f12.6)
c
c      write(*,101) ((a(i,j),i=1,n),j=1,n)
c      call invmtx2(a,n,nobs)
c      write(*,101) ((a(i,j),i=1,n),j=1,n)
c      end

       integer n, nobs,indx
       real d
       real a,a_cpy,b,b_cpy  
       dimension indx(mxlim)
       dimension a(mxlim,mxlim)         ! matrix to be inverted
       dimension a_cpy(mxlim,mxlim)     ! copy of a.  a is reused to store LU decomp
       dimension b(mxlim,mxlim)         ! Identity matrix
       dimension b_cpy(mxlim,mxlim)     ! copy of I.  b is reused in the backsub step
       
             
c  loop through making identity matrix and copy.  Also make 
c   a copy of the input matrix

      do 12 i=1,nobs
         do 11 j=1,nobs
            b(i,j)=0.0
            a_cpy(i,j)=a(i,j)
            b_cpy(i,j)=0.0

  11      continue
         b(i,i)=1.0
         b_cpy(i,i)=1.0
  12   continue

c  LU decompose the matrix.   Iterative improvement is applied twice.
c   These routines are discussed in Numerical Recipes in FORTRAN Press
c   et al pages 34-49                                                 
c   101  format(1x,5f12.6)

      call ludcmp(a,nobs,n,indx,d)       ! LU decomp
      do 13 j=1,nobs
         call lubksb(a,nobs,n,indx,b(1,j)) ! back substitute b is now Inv(a)
  13   continue

c do interative improvement 2 times.  b is the improved inverse of a.
      call mprove(a_cpy,a,nobs,n,indx,b_cpy,b)
      call mprove(a_cpy,a,nobs,n,indx,b_cpy,b)
      do 14 i=1,nobs
         do 15 j=1,nobs
         a(i,j)=b(i,j)        !return inv(a) in matrix a.
  15   continue
  14   continue
       RETURN
       end
