The variational nodal transport method is generalized for the effective treatment of multigroup criticality problems in two and three dimensions. A symbolic manipulation procedure is developed to achieve the fully automated generation of nodal response matrices in three-dimensional and non-Cartesian geometries. A red-black partitioned matrix algorithm for accelerating the solutions of the resulting within-group equations is presented, and its efficacy demonstrated. The methods are implemented as an option of the Argonne National Laboratory code DIF3D and applied to a series of five benchmark problems in x-y-z and hexagonal-z geometries. For reactors with large transport effects, the variational P3 calculations agree with accurate Monte Carlo eigenvalues to within a few hundredths to a few tenths of a percent while requiring Cray X-MP computing times ranging from tens to hundreds of seconds.