Capacitance matrix methods provide techniques for extending the use of fast Poisson solvers to arbitrary bounded regions. These techniques are further studied and developed with a focus on the three-dimensional case. A discrete analogue of classical potential theory is used as a guide in the design of rapidly convergent iterative methods. Algorithmic and programming aspects of the methods are also explored in detail. Several conjugate gradient methods are discussed for the solution of the capacitance matrix equation. A fast Poisson solver is developed which is numerically very stable even for indefinite Helmholtz equations. Variants thereof allow substantial savings in primary storage for problems on very fine meshes. Numerical results show that accurate solutions can be obtained at a cost which is proportional to that of the fast Helmholtz solver in use.