For each step of Subroutine REF , we need to solve two linear systems: (REF ) in Subroutine REF  and (REF ) in Subroutine REF . Since both of these linear systems are positive definite, we can use the preconditioned conjugate gradient (PCG) method to solve them. As (REF ) is a discretized formulation of the elliptic equation (REF ), a popular way to construct the preconditioner \(\tilde{A}_{k_m}\)  is using the multigrid (MG) method (see, e.g. [1]}), which uses the MG V-cycles associated with its coefficient matrix \(A_{k_m}\)  to approximate \(A_{k_m}^{-1}\) . For (REF ) in Subroutine REF , we still use the preconditioner \(\tilde{A}_{k_m}\)  of (REF ) for solving (REF ), though the MG V-cycles associated with \(H_{k_m}\)  may be closer to \(H_{k_m}^{-1}\) . Note that we do not use the MG V-cycles associated with \(H_{k_m}\)  because it is expensive to compute the explicit formulation of \(H_{k_m}\)  and MG V-cycles require more computation if the explicit formulation of \(H_{k_m}\)  is unknown.
