This paper presents the formulation and implementation of a high performance algorithm to compute the many-body electronic correlation energy via the random-phase approximation within density functional theory. Our approach circumvents computational inefficiencies inherent in direct approaches which exhibit quartic scaling with respect to system size. Our formulation requires solving block linear systems whose coefficient matrices are complex symmetric; these systems are of widely-varying numerical difficulty. We develop a shortterm recurrence block Krylov subspace solver for these systems and leverage a dynamic block size selection to mitigate load imbalances. This selection balances the increased cost per linear solver iteration with a reduction in the number of iterations for slowly-converging systems. Numerical experiments show that our implementation exhibits good parallel scalability, achieves faster solution times than direct approaches on even the smallest chemical system tested, and scales to larger systems and processor counts due to its cubic scaling and greater computational locality.