subroutine ludcmp(a,n,np,indx,d,nsing)
      implicit real*8 (a-h,o-z)
      parameter (nmax=1000,tiny=1.0D-15)
      dimension a(np,np),indx(n),vv(nmax)

      d = 1.0D+00
      nsing = 0

      do 12 i = 1,n
         aamax = 0.0D+00
         do 11 j = 1,n
            if (abs(a(i,j)) .gt. aamax) aamax = abs(a(i,j))
 11      continue
C         if (aamax .eq. 0.0D+00) pause 'Singular matrix'
         if (dabs(aamax) .le. 1.0D-13) then
            nsing = 1
            goto 9999
         endif
         vv(i) = 1.0D+00/aamax
 12   continue

      do 19 j = 1,n
         do 14 i = 1,j-1
            sum = a(i,j)
            do 13 k = 1,i-1
               sum = sum - a(i,k)*a(k,j)
 13         continue
            a(i,j) = sum
 14      continue
         aamax = 0.0D+00
         do 16 i = j,n
            sum = a(i,j)
            do 15 k = 1,j-1
               sum = sum - a(i,k)*a(k,j)
 15         continue
            a(i,j) = sum
            dum = vv(i)*abs(sum)
            if (dum .ge. aamax) then
               imax = i
               aamax = dum
            endif
 16      continue
         if (j .ne. imax) then
            do 17 k = 1,n
               dum = a(imax,k)
               a(imax,k) = a(j,k)
               a(j,k) = dum
 17         continue
            d = -d
            vv(imax) = vv(j)
         endif
         indx(j) = imax
         if (a(j,j) .eq. 0.0D+00) a(j,j) = tiny
         if (j .ne. n) then
            dum = 1.0D+00/a(j,j)
            do 18 i = j+1,n
               a(i,j) = a(i,j)*dum
 18         continue
         endif
 19   continue     
 
9999  return
      end

      subroutine lubksb(a,n,np,indx,b)
      implicit real*8 (a-h,o-z)
      dimension a(np,np),indx(n),b(n)

      ii = 0
 
      do 12 i = 1,n
         ll = indx(i)
         sum = b(ll)
         b(ll) = b(i)
         if (ii .ne. 0) then
            do 11 j = ii,i-1
               sum = sum - a(i,j)*b(j)
 11         continue
         else if (sum .ne. 0.0D+00) then
            ii = i
         endif
         b(i) = sum
 12   continue

      do 14 i = n,1,-1
         sum = b(i)
         if (i .lt. n) then
            do 13 j = i+1,n
               sum = sum - a(i,j)*b(j)
 13         continue
         endif
         b(i) = sum/a(i,i)
 14   continue

      return
      end

      double precision function det1(a,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=1)
      dimension a(n,n),atemp(n,n),indx(n)

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) then
         det1 = 0.0D+00
      else 
         do 11 j = 1,n
            d = d*atemp(j,j)
 11      continue
         det1 = d
      endif

      return
      end

      double precision function det3(a,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=3)
      dimension a(n,n),atemp(n,n),indx(n)

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) then
         det3 = 0.0D+00
      else 
         do 11 j = 1,n
            d = d*atemp(j,j)
 11      continue
         det3 = d
      endif

      return
      end

      double precision function det2(a,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=2)
      dimension a(n,n),atemp(n,n),indx(n)

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) then
         det2 = 0.0D+00
      else 
         do 11 j = 1,n
            d = d*atemp(j,j)
 11      continue
         det2 = d
      endif

      return
      end

      double precision function det4(a,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=4)
      dimension a(n,n),atemp(n,n),indx(n)

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) then
         det4 = 0.0D+00
      else 
         do 11 j = 1,n
            d = d*atemp(j,j)
 11      continue
         det4 = d
      endif

      return
      end

      double precision function det8(a,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=8)
      dimension a(n,n),atemp(n,n),indx(n)

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) then
         det8 = 0.0D+00
      else 
         do 11 j = 1,n
            d = d*atemp(j,j)
 11      continue
         det8 = d
      endif

      return
      end

      subroutine inv8(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=8)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv200(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=200)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv105(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=105)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv36(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=36)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv23(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=23)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end
                 
      subroutine inv22(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=22)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv29(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=29)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv27(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=27)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv54(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=54)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv51(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=51)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv28(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=28)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv30(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=30)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv15(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=15)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv19(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=19)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv3(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=3)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv2(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=2)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv1(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=1)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv11(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=11)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv16(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=16)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv20(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=20)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv6(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=6)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv4(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=4)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv49(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=49)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv97(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=97)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv5(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=5)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv7(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=7)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv10(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=10)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv9(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=9)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv13(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=13)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv17(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=17)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv12(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=12)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv14(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=14)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv18(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=18)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv21(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=21)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine inv37(a,b,nsing)
      implicit real*8 (a-h,o-z)
      parameter (n=37)
      dimension a(n,n),b(n,n),atemp(n,n),eye(n,n),indx(n)

      call inimat(eye,n,n)
      do 10 i = 1,n
         eye(i,i) = 1.0D+00
 10   continue

      call copmat(a,atemp,n,n)

      call ludcmp(atemp,n,n,indx,d,nsing)
      if (nsing .eq. 1) goto 9999
      do 11 j = 1,n
         call lubksb(atemp,n,n,indx,eye(1,j))
 11   continue

      call copmat(eye,b,n,n)

9999  return
      end

      subroutine scalar(a,b,c,n,k)
      implicit real*8 (a-h,o-z)
      dimension a(n,k),b(n,k)
      
      do 10 i = 1,n
         do 20 j = 1,k
            b(i,j) = c*a(i,j)
 20      continue
 10   continue

      return
      end

      subroutine inimat(a,n,k)
      implicit real*8 (a-h,o-z)
      dimension a(n,k)

      do 10 i = 1,n
         do 20 j = 1,k
            a(i,j) = 0.0D+00
 20      continue
 10   continue
  
      return
      end
 
      subroutine copmat(a,b,n,k)
      implicit real*8 (a-h,o-z)
      dimension a(n,k),b(n,k)

      do 10 i = 1,n
         do 20 j = 1,k
            b(i,j) = a(i,j)
 20      continue
 10   continue

      return
      end

      subroutine add2(a,b,c,n,k)
      implicit real*8 (a-h,o-z)
      dimension a(n,k),b(n,k),c(n,k)

      do 10 i = 1,n
         do 20 j = 1,k
            c(i,j) = a(i,j) + b(i,j)
 20      continue
 10   continue
 
      return
      end

      subroutine sub2(a,b,c,n,k)
      implicit real*8 (a-h,o-z)
      dimension a(n,k),b(n,k),c(n,k)

      do 10 i = 1,n
         do 20 j = 1,k
            c(i,j) = a(i,j) - b(i,j)
 20      continue
 10   continue
 
      return
      end

      subroutine add3(a,b,c,d,n,k,temp)
      implicit real*8 (a-h,o-z)
      dimension a(n,k),b(n,k),c(n,k),d(n,k),temp(n,k)

      call add2(a,b,temp,n,k)
      call add2(c,temp,d,n,k)

      return
      end

      subroutine tr(a,b,n,k)
      implicit real*8 (a-h,o-z)
      dimension a(n,k),b(k,n)

      do 10 i = 1,n
         do 20 j = 1,k
            b(j,i) = a(i,j)
 20      continue
 10   continue

      return
      end

C      subroutine mult2(a,b,c,n,k,l)
C      implicit real*8 (a-h,o-z)
C      dimension a(n,k),b(k,l),c(n,l)
C
C      do 10 i = 1,n
C         do 11 j = 1,l
C            c(i,j) = 0.0D+00
C            do 12 m = 1,k
C               c(i,j) = c(i,j) + a(i,m)*b(m,j)
C 12         continue
C 11      continue
C 10   continue
C
C      return
C      end

      subroutine mult2(a,b,c,n,k,l)
      implicit real*8 (a-h,o-z)
      dimension a(n,k),b(k,l),c(n,l)

      if (l .eq. 1) then 
         do 20 i = 1,n
            c(i,1) = 0.0D+00
            do 22 m = 1,k
               c(i,1) = c(i,1) + a(i,m)*b(m,1)
 22         continue
 20      continue
      else
         do 10 i = 1,n
            do 11 j = 1,l
               c(i,j) = 0.0D+00
               do 12 m = 1,k
                  c(i,j) = c(i,j) + a(i,m)*b(m,j)
 12            continue
 11         continue
 10      continue
      endif

      return
      end

      subroutine mult3(a,b,c,d,n,k,l,m,temp)
      implicit real*8 (a-h,o-z)
      dimension a(n,k),b(k,l),c(l,m),d(n,m),temp(n,l)

      call mult2(a,b,temp,n,k,l)
      call mult2(temp,c,d,n,l,m)

      return
      end

      subroutine dkron(x,n1,m1,y,n2,m2,z)
	implicit real*8(a-h,o-z)
      dimension x(n1,m1),y(n2,m2),z(n1*n2,m1*m2)

	do 13 i = 1,m1
	   do 12 j = 1,n1
	      do 11 ii = 1,m2
	         do 10 jj = 1,n2
	            z(n1*(i-1)+ii,m1*(j-1)+jj) = x(i,j)*y(ii,jj)
 10            continue
 11         continue
 12      continue
 13   continue

      return
	end