A multilevel accelerated lambda iteration (MALI) method for radiative transfer calculations with partial frequency redistribution (PRD) is presented. The method, which is based on Rybicki & Hummer's complete frequency redistribution (CRD) formalism with full preconditioning, consistently accounts for overlapping radiative transitions. Its extension to PRD is implemented in a very natural way through the use of the Ψ operator operating on the emissivity rather than the commonly used Λ operator, which operates on the source function. Apart from requiring an additional inner computational loop to evaluate the PRD emission-line profiles with fixed population numbers, implementation of the presented method requires only a trivial addition of computer code. Since the presented method employs a diagonal operator, it is easily extended to different geometries. Currently, it has been implemented for one-, two-, and three-dimensional Cartesian grids and spherical symmetry. In all cases, the speed of convergence with PRD is very similar to that in CRD, with the former sometimes even surpassing the latter. Sample calculations exhibiting the favorable convergence behavior of the PRD code are presented in the case of the Ca II H and K lines, the Mg II h and k lines, and the hydrogen Lyα and Lyβ lines in a one-dimensional solar model and the Ca II resonance lines in a two-dimensional flux-sheet model.