Consider the NxN system of linear equations M x = b, where the coefficient matrix M is large, sparse, and nonsymmetric. Assume that M can be factored in the form M = L D U, where L is a lower triangular matrix, D is a diagonal matrix, and U is a unit upper triangular matrix. Such systems arise frequently in scientific computation, e.g., in finite difference and finite element approximations to non-self-adjoint elliptic boundary value problems. This report presents a package of efficient, reliable, well-documented, and portable FORTRAN subroutines for solving these systems.