A numerical method for the solution of low Reynolds number flow (Stokes flow) past particles of arbitrary shape is described. The method uses a predictor–corrector approach, with a rough initial calculation using the Dabros method followed by an application of the completed double layer boundary integral equation method (CDL‐BIEM). The method is especially well suited for particles with very rough and uneven surfaces. The mobility and stresslet of a sphere with a spiked surface are computed, to illustrate the utility of the new method. For spheres with small protrusions, our results show a small change in the stresslet and particle mobility. For protrusions on the order of the sphere radius, there is a significant increase in the stresslet, leading to Einstein viscosity coefficients on the order of 5, compared to 5/2 for the perfect sphere. We also examine algorithms and the underlying theory behind conjugate gradient and conjugate residual methods that appear to work well with the larger CDL‐BIEM problems. For simpler geometries, the discretized equations may be readily handled by a number of equation solvers, including direct methods such as LU factorization. For large (50 000 by 50 000 and greater) dense linear systems associated with multiparticle problems in complex geometries, the iterative methods show great promise, especially on high‐performance parallel computer architectures.