Descriptions of implicit solvent models like the Poisson-Boltzmann equation express the electrostatic potential as the solution of an elliptic PDE with delta function source terms. These equations are particularly good fits with iterative solvers that use a fast N-body solver as a preconditioner. Also, modern architectures favor such algorithms due to their high ratio of floating-point operations to memory references. Two types of N-body solvers can be distinguished: hierarchical-clustering algorithms, such as the celebrated fast multipole method, and kernel-splitting algorithms, such as the popular particle--mesh Ewald method. By formulating the problem as that of computing a matrix--vector product, the basic structure of these algorithms is elucidated. Additionally, evidence is presented indicating that kernel-splitting algorithms are much to be preferred for molecular simulations and that the virtually unknown multilevel summation method of Brandt and Lubrecht is the best among these. This method uses hierarchical interpolation of interaction potentials on nested grids to calculate energies and forces in linear time for both periodic and nonperiodic boundary conditions. This is joint work with David Hardy. Its application to the Poisson(-Boltzmann) equation is being pursued in a collaboration with Stephen Bond.