  A fast direct inversion scheme for the large sparse systems of linear
equations resulting from the discretization of elliptic partial differential
equations in two dimensions is given. The scheme is described for the
particular case of a discretization on a uniform square grid, but can be
generalized to more general geometries. For a grid containing $N$ points, the
scheme requires $O(N \log^{2}N)$ arithmetic operations and $O(N \log N)$
storage to compute an approximate inverse. If only a single solve is required,
then the scheme requires only $O(\sqrt{N} \log N)$ storage; the same storage is
sufficient for computing the Dirichlet-to-Neumann operator as well as other
boundary-to-boundary operators. The scheme is illustrated with several
numerical examples. For instance, a matrix of size $10^6 \times 10^6$ is
inverted to seven digits accuracy in four minutes on a 2.8GHz P4 desktop PC.
