  Integral equation methods for solving the Laplace-Beltrami equation on the
unit sphere in the presence of multiple "islands" are presented. The surface of
the sphere is first mapped to a multiply-connected region in the complex plane
via a stereographic projection. After discretizing the integral equation, the
resulting dense linear system is solved iteratively using the fast multipole
method for the 2D Coulomb potential in order to calculate the matrix-vector
products. This numerical scheme requires only O(N) operations, where $N$ is the
number of nodes in the discretization of the boundary. The performance of the
method is demonstrated on several examples.
