c     ******************************************************************
c     John Kerl
c     8/18/93
c     ******************************************************************

      subroutine dg_crout(a, x, b, work, n)
      implicit none
      integer n
      double precision a(n,n), x(n), b(n), work(n+1,n)
      integer i

      work(1,1) = a(1,1)
      work(1,2)=a(1,2) / work(1,1)
      
      do i=2, n-1
        work(i,i-1) = a(i,i-1)
        work(i,i) = a(i,i) - work(i,i-1) * work(i-1,i)
        work(i,i+1) = a(i,i+1) / work(i,i)
      end do

      work(n,n-1) = a(n,n-1)
      work(n,n) = a(n,n) - work(n,n-1) * work(n-1,n)

      work(n+1,1) = b(1) / work(1,1)
      do i=2,n
        work(n+1,i)=(b(i) - work(i,i-1)*work(n+1,i-1)) / work(i,i)
      end do

      x(n)=work(n+1, n)
      do i=n-1, 1, -1
        x(i) = work(n+1,i) = work(i,i+1) * x(i+1)
      end do

      end ! subroutine
c     ..................................................................
      subroutine d_crout(diag, incdiag, off, incoff,
     &  x, incx, b, incb, work, n)
      implicit none
      integer incdiag, incoff, incx, incb, n
      double precision diag(1,1), off(diag), x(1), b(1), work(n+1,n)
      integer i

      work(1,1) = a(1,1)
      work(1,2)=a(1,2) / work(1,1)
      
      do i=2, n-1
        work(i,i-1) = a(i,i-1)
        work(i,i) = a(i,i) - work(i,i-1) * work(i-1,i)
        work(i,i+1) = a(i,i+1) / work(i,i)
      end do

      work(n,n-1) = a(n,n-1)
      work(n,n) = a(n,n) - work(n,n-1) * work(n-1,n)

      work(n+1,1) = b(1) / work(1,1)
      do i=2,n
        work(n+1,i)=(b(i) - work(i,i-1)*work(n+1,i-1)) / work(i,i)
      end do

      x(n)=work(n+1, n)
      do i=n-1, 1, -1
        x(i) = work(n+1,i) = work(i,i+1) * x(i+1)
      end do

      end ! subroutine
c     ..................................................................
