Least-squares methodology applied to the 3-D neutron transport equation is a practical alternative to current numerical methods. The least-squares transport method, through a properly scaled least-squares functional, defines a minimization principle which allows easy construction of the solution procedure. The scaling is chosen to facilitate an approximation of the correct behavior of the solution in the different parameter regimes.Representing the angular dependence with a finite expansion of spherical harmonics creates a self-adjoint, positive definite system of moment equations to be spatially discretized with finite element functions. To approximately invert this system of equations, a conjugate gradient method with a block Jacobi preconditioner, that uses multigrid to invert the diagonal blocks, is defined and discussed. For nearly all parameter regimes and all diagonal blocks, standard multigrid is sufficient as a preconditioner. However, for certain parameter regimes, multigrid performs poorly on the diagonal blocks associated with the first-order moments. A new multigrid algorithm is proposed which simultaneously addresses the three first-order moments and recovers essential multilevel properties.