  For analyzing maps of the cosmic microwave background sky, it is necessary to
mask out the region around the galactic equator where the parasitic foreground
emission is strongest as well as the brightest compact sources. Since many of
the analyses of the data, particularly those searching for non-Gaussianity of a
primordial origin, are most straightforwardly carried out on full-sky maps, it
is of great interest to develop efficient algorithms for filling in the missing
information in a plausible way. We explore practical algorithms for filling in
based on constrained Gaussian realizations. Although carrying out such
realizations is in principle straightforward, for finely pixelized maps as will
be required for the Planck analysis a direct brute force method is not
numerically tractable. We present some concrete solutions to this problem, both
on a spatially flat sky with periodic boundary conditions and on the pixelized
sphere. One approach is to solve the linear system with an appropriately
preconditioned conjugate gradient method. While this approach was successfully
implemented on a rectangular domain with periodic boundary conditions and
worked even for very wide masked regions, we found that the method failed on
the pixelized sphere for reasons that we explain here. We present an approach
that works for full-sky pixelized maps on the sphere involving a kernel-based
multi-resolution Laplace solver followed by a series of conjugate gradient
corrections near the boundary of the mask.
