  A novel and scalable geometric multi-level algorithm is presented for the
numerical solution of elliptic partial differential equations, specially
designed to run with high occupancy of streaming processors inside Graphics
Processing Units(GPUs). The algorithm consists of iterative, superposed
operations on a single grid, and it is composed of two simple full-grid
routines: a restriction and a coarsened interpolation-relaxation. The
restriction is used to collect sources using recursive coarsened averages, and
the interpolation-relaxation simultaneously applies coarsened finite-difference
operators and interpolations. The routines are scheduled in a saw-like refining
cycle. Convergence to machine precision is achieved repeating the full cycle
using accumulated residuals and successively collecting the solution. Its total
number of operations scale linearly with the number of nodes. It provides an
attractive fast solver for Boundary Value Problems (BVPs), specially for
simulations running entirely in the GPU. Applications shown in this work
include the deformation of two-dimensional grids, the computation of
three-dimensional streamlines for a singular trifoil-knot vortex and the
calculation of three-dimensional electric potentials in heterogeneous
dielectric media.
