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.