First-order least-squares method of a distributed optimal control problem for the incompressible Stokes equations is considered. An optimality system for the optimal solution are reformulated to the equivalent first-order system by introducing the vorticity and then the least-squares functional corresponding to the system is defined in terms of the sum of the squared $H^{-1}$ and $L^2$ norms of the residual equations of the system. Finite element approximations are studied and optimal error estimates are obtained. Resulting linear system of the optimality system is symmetric and positive definite. The V-cycle multigrid method is applied to the system to test computational efficiency.