Multigrid-preconditioned Krylov methods are applied to within-group response matrix equations of the type derived from the variational nodal method for neutron transport with interface conditions represented by orthogonal polynomials in space and spherical harmonics in angle. Since response matrix equations result in nonsymmetric coefficient matrices, the generalized minimal residual (GMRES) Krylov method is employed. Two acceleration methods are employed: response matrix aggregation and multigrid preconditioning. Without approximation, response matrix aggregation combines fine-mesh response matrices into coarse-mesh response matrices with piecewise-orthogonal polynomial interface conditions; this may also be viewed as a form of nonoverlapping domain decomposition on the coarse grid. Two-level multigrid preconditioning is also applied to the GMRES method by performing auxiliary iterations with one degree of freedom per interface that conserve neutron balance for three types of interface conditions: (a) p preconditioning is applied to orthogonal polynomial interface conditions (in conjunction with matrix aggregation), (b) h preconditioning to piecewise-constant interface conditions, and (c) h-p preconditioning to piecewise-orthogonal polynomial interface conditions. Alternately, aggregation is employed outside the GMRES algorithm to coarsen the grid, and multigrid preconditioning is then applied to the coarsened equations. The effectiveness of the combined aggregation and preconditioning techniques is demonstrated in two dimensions on a fixed-source, within-group neutron diffusion problem approximating the fast group of a pressurized water reactor configuration containing six fuel assemblies.