The evaluation of a domain integral is the dominant bottleneck in the numerical solution of viscous flow problems by vorticity methods, which otherwise demonstrate distinct advantages over primitive variable methods. By applying a Bames-Hut multipole acceleration technique, the operation count for the integration is reduced from O(N-2) to O(N log N), while the memory requirements are reduced from O(N-2) to O(N). The algorithmic parameters that are necessary to achieve such scaling are described. The parallelization of the algorithm is crucial if the method is to be applied to realistic problems. A parallelization procedure which achieves almost perfect scaling is shown. Finally, numerical experiments on a driven cavity benchmark problem are performed. The actual increase in performance and reduction in storage requirements match theoretical predictions well, and the scalability of the procedure is very good. Copyright (C) 2003 John Wiley Sons, Ltd.