In this paper a boundary element method based on primitive variable formulation for two-dimensional steady Stokes flow is presented. By the use of Lorentz, reciprocal theorem for steady incompressible viscous flow the boundary integral equations for the velocity field are derived. Discretization of boundary of the problem leads to a system of equations wherein all the unknown velocity and traction components are on the boundary nodes. Once these boundary velocities and tractions are calculated, the velocities, pressure, shearing stresses and vorticity at each internal point can be obtained from the boundary integral representations along with various kernel functions. Two examples including flow squeezed between two parallel plates as well as wall driven cavity flow are employed to test the algorithms. It is found that excellent agreement is obtained between the boundary element results and the available analytical and finite element results.