A fast iterative algorithm is presented for the numerical solution of large linear systems that are encountered in multiparticle Stokes flows. It is applicable to solid particles of arbitrary shape, and finds the translation and rotation velocities when total forces and torques acting on these particles are given (mobility problems). An exact result for the stresslet is also given. The method is based on recently developed boundary integral equations, “the canonical equations for mobility and resistance problems.” These are well‐posed Fredholm equations of the second kind, for mobility problems or problems with arbitrary velocity boundary conditions. They are modified for the direct iterative solution of mobility problems, leading to fast numerical computations. For a single sphere the iteration operator is spectrally decomposed analytically. The convergence rate of the iterations is deduced, and supporting numerical observations are presented. Fast rate of convergence is numerically observed for multisphere problems even at small separations. Simulations in polymeric and suspension rheology and protein chemistry are envisioned as potential applications. Numerous solutions are needed when simulating such systems over some period of time. With nearly equal consecutive boundary configurations, each previous solution provides a good initial guess for the next iterative solution.