A common challenge in simulating dense suspension of rigid particles in Stokes flow is the numerical inaccuracies and instabilities that arises due to particle collisions. To overcome this problem, a strong repulsive potential between particles is often prescribed. This in turn leads to numerical stiffness and dramatic reduction in stable time-step sizes. In this work, we eliminate such stiffness by introducing contact constraints explicitly and solving the hydrodynamic equations in tandem with a linear complementarity problem with inequality constraints. Satisfaction of Newtons third law for the collision force is explicitly guaranteed, allowing the consistent calculation of collision stresses. Efficient parallelization for shared-memory and distributed-memory architectures is also implemented. This method can be coupled to any Stokes hydrodynamics solver for particles with various shapes and allows us to simulate $10^4-10^7$ spheres on a laptop, depending on the cost of the Stokes hydrodynamics solver. We demonstrate its performance on a range of applications of rigid suspensions.