An iterative algorithm based on fast Fourier transforms is proposed to solve the Poisson equation for systems of heterogeneous permittivity (e.g., solute cavity in a solvent) under periodic boundary conditions. The method makes explicit use of the dipole-dipole interaction tensor, and is thus easily generalizable to arbitrary forms of electrostatic interactions (e.g., Coulomb's law with straight or smooth cutoff truncation). The convergence properties of the algorithm and the influence of various model parameters are investigated in detail, and a set of appropriate values for these parameters is determined. The algorithm is further tested by application to three types of systems (a single spherical ion, two spherical ions, and small biomolecules), and comparison with analytical results (single ion) and with results obtained using a finite-difference solver under periodic boundary conditions. The proposed algorithm performs very well in terms of accuracy and convergence properties, with an overall speed comparable in the current implementation to that of a typical finite-difference solver. Future developments and applications of the algorithm will include: (i) the assessment of periodicity- and cutoff-induced artifacts in explicit-solvent simulations; (ii) the design of new electrostatic schemes for explicit-solvent simulations mimicking more accurately bulk solution; (iii) a faster evaluation of solvation free energies based on continuum electrostatics in cases where periodicity-induced artifacts can be neglected.