Direct and indirect boundary element methods, accelerated via the fast multipole method, are applied to numerical simulation of room acoustics for rooms of volume ∼150 m3 and frequencies up to 5 kHz on a workstation. As the parameter kD (wavenumber times room diameter) is large, stabilization of the previously developed fast algorithms is required for accuracy. A stabilization scheme is one of the key contributions of this paper. The computations are validated using well-known image source solutions for shoebox shaped rooms. Computations for L-shaped rooms are performed to illustrate the ability to capture diffractions. The ability to model in-room baffles and boundary openings (doors/windows) is also demonstrated. The largest case has kD > 1100 with a discretization of size 6 × 106 elements. The performance of different boundary integral formulations was compared, and their rates of convergence using a preconditioned flexible generalized minimal residual method were found to be substantially different. These promising results suggest a path to efficient simulations of room acoustics via high performance boundary element methods.