  We present an iterative scheme, reminiscent of the Multigrid method, to solve
large boundary value problems with Probabilistic Domain Decomposition (PDD). In
it, increasingly accurate approximations to the solution are used as control
variates in order to reduce the Monte Carlo error of the following
iterates--resulting in an overall acceleration of PDD for a given error
tolerance. The key ingredient of the proposed algorithm is the ability to
approximately predict the speedup with little computational overhead and in
parallel. Besides, the theoretical framework allows to explore other aspects of
PDD, such as stability. One numerical example is worked out, yielding an
improvement of between one and two orders of magnitude over the previous
version of PDD.
