  We present a fast direct solver for structured linear systems based on
multilevel matrix compression. Using the recently developed interpolative
decomposition of a low-rank matrix in a recursive manner, we embed an
approximation of the original matrix into a larger, but highly structured
sparse one that allows fast factorization and application of the inverse. The
algorithm extends the Martinsson/Rokhlin method developed for 2D boundary
integral equations and proceeds in two phases: a precomputation phase,
consisting of matrix compression and factorization, followed by a solution
phase to apply the matrix inverse. For boundary integral equations which are
not too oscillatory, e.g., based on the Green's functions for the Laplace or
low-frequency Helmholtz equations, both phases typically have complexity O(N)
in two dimensions, where $N$ is the number of discretization points. In our
current implementation, the corresponding costs in three dimensions are
$O(N^{3/2})$ and $O(N \log N)$ for precomputation and solution, respectively.
Extensive numerical experiments show a speedup of $\sim 100$ for the solution
phase over modern fast multipole methods; however, the cost of precomputation
remains high. Thus, the solver is particularly suited to problems where large
numbers of iterations would be required. Such is the case with ill-conditioned
linear systems or when the same system is to be solved with multiple right-hand
sides. Our algorithm is implemented in Fortran and freely available.
