subroutine decomp ( n, ndim, a, ip )

c*********************************************************************72
c
c  matrix triangularization by gaussian elimination.
c
c  input..
c
c    n = order of matrix.
c    ndim = declared dimension of array a.
c    a = matrix to be triangularized.
c
c  output..
c
c    a(i,j), i.le.j = upper triangular factor u.
c    a(i,j), i.gt.j = multipliers = lower triangular
c                     factor, i-l.
c    ip(k), k.lt.n = index of k-th pivot row.
c    ip(n) = (-1)**(number of interchanges) or 0.
c
c  use 'solve' to obtain solution of linear system.
c  determ(a) = ip(n)*a(1,1)*a(2,2)*...*a(n,n).
c  if ip(n) = 0, a is singular, solve will divide by zero.
c  interchanges finished in u, only partly in l.
c
      real a(ndim,ndim), t
      integer ip(ndim)

      ip(n) = 1

      do k = 1, n

        if ( k .eq. n ) go to 5
        kp1 = k + 1
        m = k
        do i = k + 1, n
          if ( abs ( a(i,k)) .gt. abs ( a(m,k) ) ) m = i
        end do
        ip(k) = m
        if ( m .ne. k ) ip(n) = -ip(n)
        t = a(m,k)
        a(m,k) = a(k,k)
        a(k,k) = t
        if ( t .eq. 0. ) go to 5
        do i = k + 1, n
          a(i,k) = - a(i,k) / t
        end do

        do j = k + 1, n
          t = a(m,j)
          a(m,j) = a(k,j)
          a(k,j) = t
          do i = k + 1, n
            a(i,j) = a(i,j) + a(i,k) * t
          end do
        end do

5       if ( a(k,k) .eq. 0.0 ) ip(n) = 0

      end do

      return
      end
      subroutine solve ( n, ndim, a, b, ip )

c*********************************************************************72
c
cc solve() solves a linear system, a*x = b.
c
c  input..
c
c    n = order of matrix.
c    ndim = declared dimension of array a.
c    a = triangularized matrix obtained from 'decomp'.
c    b = right hand side vector.
c    ip = pivot vector obtained from 'decomp'.
c  do not use if decomp has set ip(n) = 0.
c
c  output..
c
c    b = solution vector, x.
c
      real a(ndim,ndim)
      real b(ndim)
      integer ip(ndim)
      real t

      do k = 1, n - 1
        kp1 = k + 1
        m = ip(k)
        t = b(m)
        b(m) = b(k)
        b(k) = t
        do i = k + 1, n
          b(i) = b(i) + a(i,k) * t
        end do
      end do

      do kb = 1, n - 1
        km1 = n - kb
        k = km1 + 1
        b(k) = b(k) / a(k,k)
        t = - b(k)
        do i = 1, km1
          b(i) = b(i) + a(i,k) * t
        end do
      end do

      b(1) = b(1) / a(1,1)

      return
      end