We present a Multilevel Quasi-Monte Carlo algorithm for the solution of elliptic partial differential equations with random coefficients and inputs. By combining the multilevel sampling idea with randomly shifted rank-1 lattice rules, the algorithm constructs an estimator for the expected value of some functional of the solution. The error analysis of this estimator provides a formula for the optimal number of samples at each level. The efficiency of the method is illustrated on a three-dimensional subsurface flow problem with lognormal diffusion coefficient. This example is particularly challenging because of the small correlation length, and thus the large number of uncertainties that must be included. We numerically show that it is possible to achieve a cost almost inversely proportional to the requested tolerance on the root-mean-square error.