A multigrid method for the solution of finite difference approximations of elliptic PDEs is introduced. A parallelizable version of it, suitable for two and multi level analysis, is also defined, and serves as a theoretical tool for deriving an optimal implementation for the main version. For indefinite Helmholtz equations, this analysis provides a prediction of the optimal mesh size for the coarsest grid used. Numerical experiments show the applicability of the method to 3-d diffusion problems with discontinuous coefficients and highly indefinite Helmholtz equations.