A fast direct solver for network matrices

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.