      program lutest
      implicit none
      integer dim, n
      parameter (dim=4)
      double precision  a(dim,dim), l(dim,dim), 
     &     u(dim,dim), x(dim), b(dim), z(dim)
      integer ier

      read (5, *) n
      if (n .gt. dim) then
         write (0,*) "Matrix dimensions too large;"
         write (0,*) "  please adjust and re-compile lutest.f."
         stop
      end if
      call d_read_matrix(a, n, n, ier)
      call d_read_matrix(b, 1, n, ier)
      if (ier .ne. 0) call usage

      call d_lusolve(a, x, b, z, l, u, n, ier)
      if ( ier .ne. 0) then
         write (0,*) "lutest:  Couldn't solve system of equations."
      else
         call dg_write_matrix(x, n, 1, '', 0)
      end if

      end ! program lutest

      subroutine usage()
      write (0,*) "lutest:  Input is n (dimension of square system),"
      write (0,*) "   A (matrix), b (right-hand-side vector)."
      return
      end ! subroutine usage
