We present a new least-squares finite element method for the steady Oldroyd type viscoelastic fluids. The overall iterative procedure combines a nonlinear nested iteration, where adaptive mesh refinement is based on a nonlinear least-squares functional. Each linear step is solved by a least-squares finite element minimization. The homogeneous least-squares functional is shown to be equivalent to a natural norm, and, under sufficient smoothness assumptions, finite element error bounds are shown to be optimal when using conforming piecewise polynomial elements for the velocity, pressure and extra stress, and Raviart-Thomas finite elements for the total stress. In the absence of full regularity, a local weighted-norm approach is used to remove effects of corner singularities. Numerical results are given for an Oldroyd-B fluid in a 4:1 contraction, showing optimal reduction of the least-squares functional. (C) 2009 Elsevier B.V. All rights reserved.