Clusters containing bubbles and solid particles are used in many fields of industry. For instance, the study of bubble-particle interaction can be useful for surface cleaning in microelectronics and froth flotation in oil distillation. The present work discusses the joint 3D dynamics of bubbles and solid spherical particles in the presence of an acoustic field in an unbounded ideal incompressible liquid. To solve the problem, we chose the boundary element method (BEM) which requires only the discretization of the boundary of the computational domain. However, the application of the conventional BEM for the direct simulation of large particle-bubble systems is normally limited by memory, computational complexity, and speed. To perform such simulations, we propose a numerical approach based on the GPU acceleration of the BEM code. The method includes the solution of linear algebraic equations with a dense matrix of special type, using for this the generalized minimal residual method with calculation of the matrix-vector product on graphics processors using CUDA technology. We also discuss the performance of the method and run test computations of bubble-particle clusters.