Pseudo-code for the MG V-cycle. The matrices involved are defined in the main text. Basic-V-Cycle: the clean textbook version, exposing the basic structure of the algorithm. The important feature of the algorithm is that the solution vector x is never transferred directly between levels. Instead, a residual r_H is computed, which takes the role as the right-hand side b on the coarser level. The coarse solution is a correction c_H that is then added to the solution vector x. The resulting full V-cycle is a symmetric linear operator that can be used as a preconditioner for CG. We keep recursing until there are less than 1000 coefficients left, and then solve using a pseudo-inverse based on the SVD, , as G_h may, depending on the size of the mask, be singular due to D_h being truncated at L_h. Optimized-V-Cycle: in this code we have inserted and , and then reorganized the expressions so that the restriction and interpolation steps each share one SHT with the corresponding application of G_h. The J operator denotes . The effect of this operator is to zero out any contribution that falls outside of the mask in the (full sky) spherical harmonic vectors; but numerical experiments indicate that the term can in practice be neglected also when using a small mask. In our numerical experiments we approximate J ≈ I, reducing the total number of SHTs to six per level. The comments indicate the required resolution for each SHT, with L_h (L_H) referring to fine (coarse) harmonic band-limit, and N_h (N_H) referring to fine (coarse) grid.