Multigrid method and acceleration by conjugate gradient method for first-order system least squares (FOSLS) using bilinear finite elements are developed for various boundary value problems of planar linear elasticity. They are two-stage algorithms that first solve for the displacement flux variable, then for the displacement itself. This paper focuses on solving for the displacement flux variable only. Numerical results show that the convergence is uniform even as the material becomes nearly incompressible. Computations for convergence factors and discretization errors are included. Heuristic arguments to improve the convergences are discussed as well.