We develop a rigid multiblob method for numerically solving the mobility problem for suspensions ofpassive and active rigid particles of complex shape in Stokes flow in unconfined, partially confined,and fully confined geometries. As in a number of existing methods, we discretize rigid bodies using acollection of minimally resolved spherical blobs constrained to move as a rigid body, to arrive at apotentially large linear system of equations for the unknown Lagrange multipliers and rigid-bodymotions. Here we develop a block-diagonal preconditioner for this linear system and show that astandard Krylov solver converges in a modest number of iterations that is essentially independentof the number of particles. Key to the efficiency of the method is a technique for fastcomputation of the product of the blob-blob mobility matrix and a vector. For unboundedsuspensions, we rely on existing analytical expressions for the Rotne–Prager–Yamakawatensor combined with a fast multipole method (FMM) to obtain linear scaling in thenumber of particles. For suspensions sedimented against a single no-slip boundary, weuse a direct summation on a graphical processing unit (GPU), which gives quadraticasymptotic scaling with the number of particles. For fully confined domains, such as periodicsuspensions or suspensions confined in slit and square channels, we extend a recentlydeveloped rigid-body immersed boundary method by B. Kallemov, A. P. S. Bhalla,B. E. Griffith, and A. Donev (Commun. Appl. Math. Comput. Sci.11 (2016), no. 1,79–141) to suspensions of freely moving passive or active rigid particles at zero Reynoldsnumber. We demonstrate that the iterative solver for the coupled fluid and rigid-bodyequations converges in a bounded number of iterations regardless of the system size. In ourapproach, each iteration only requires a few cycles of a geometric multigrid solver for thePoisson equation, and an application of the block-diagonal preconditioner, leading to linearscaling with the number of particles. We optimize a number of parameters in the iterativesolvers and apply our method to a variety of benchmark problems to carefully assess theaccuracy of the rigid multiblob approach as a function of the resolution. We also model thedynamics of colloidal particles studied in recent experiments, such as passive boomerangs in aslit channel, as well as a pair of non-Brownian active nanorods sedimented against awall.