      subroutine d_lusolve(A, x, b, z, L, U, n, ier)
      implicit none
      integer n, ier
      double precision A(n,n), x(n), b(n), L(n,n), U(n,n), z(n)
      logical debug
      debug=.true.

      call d_ludecomp(A, L, U, n, ier)
      if (ier .eq. 1) then
         write(0, *) "d_lusolve: Matrix couldn't be factored."
         ier=1
         return
      end if
      if (debug) then
        write (*,*) 'A is:'
        call d_write_matrix(A, n, n, '', 0, 'g', -1)
        write (*,*) 'L is:'
        call d_write_matrix(L, n, n, '', 0, 'lt', -1)
        write (*,*) 'U is:'
        call d_write_matrix(U, n, n, '', 0, 'ut', -1)
      end if
        
      call d_backsubplus(l, b, z, n, ier)
      if (ier .eq. 1) then
         write(0, *)
     &"d_lusolve: Could not back-substitute solution (division by zero)."
         ier=1
         return
      end if
      call d_backsubminus(u, z, x, n, ier)
      if (ier .eq. 1) then
         write(0, *) 
     &"d_lusolve: Could not back-substitute solution (division by zero)."
         ier=1
         return
      end if

      ier=0
      end ! subroutine d_lusolve

c     ******************************************************************
      subroutine d_ludecomp(A, L, U, n, ier) 
                                ! Need to pack L,U into one mat
      implicit none
      integer n, ier
      double precision A(n,n), L(n,n), U(n,n)
      integer i, j, k
      
c     Initialize workspace to the identity matrix
      do i=1,n
         do j=1,n
            L(i,j)=0.0
            U(i,j)=0.0
         end do
      end do
      do i=1,n
         L(i,i)=1.0
      end do
      
      U(1,1)=A(1,1)
      if (U(1,1)*L(1,1) .eq. 0.0) then
         write (0,*) "d_ludecomp: Factorization impossible."
         ier=1
         return
      end if
      do j=2,n
         U(1,j) = A(1,j) / L(1,1)
         L(j,1) = A(j,1) / U(1,1)
      end do

      do i=2, n-1
         U(i,i) = A(i,i)
         do k=1, i-1
            U(i,i) = U(i,i) - L(i,k) * U(k,i)
         end do  ! k=1, i-1
         if (L(i,i)*U(i,i) .eq. 0.0) then
            write (0,*) "d_ludecomp: Factorization impossible"
            ier = 1
            return
         end if
         do j=i+1, n
            U(i,j) = A(i,j)
            do k=1, i-1
               U(i,j) = U(i,j) - L(i,k) * U(k,j)
            end do 
            U(i,j) = U(i,j) / L(i,i)
            L(j,i) = A(j,i)
            do k=1, i-1
               L(j,i) = L(j,i) - L(j,k) * U(k,i)
            end do
            L(j,i) = L(j,i) / U(i,i)
         end do ! j=i+1, n
      end do ! i=2, n-1

      U(n,n) = A(n,n)
      do k=1, n-1
         U(n,n) = U(n,n) - L(n,k)*U(k,n)
      end do

      ier=0
      end ! subroutine d_ludecomp
