Popular eigensolvers such as block-Lanczos require repeated inversion of an eigenmatrix. This is a bottleneck in large-scale modal problems with millions of degrees of freedom. On the other hand, the classic Rayleigh–Ritz conjugate gradient method only requires a matrix-vector multiplication, and is therefore potentially scalable to such problems. However, as is well-known, the Rayleigh–Ritz has serious numerical deficiencies, and has largely been abandoned by the finite-element community. In this paper, we address these deficiencies through subspace augmentation, and consider a subspace augmented Rayleigh–Ritz conjugate gradient method (SaRCG). SaRCG is numerically stable and does not entail explicit inversion. As a specific application, we consider the modal analysis of geometrically complex structures discretized via nonconforming voxels. The resulting large-scale eigenproblems are then solved via SaRCG. The voxelization structure is also exploited to render the underlying matrix-vector multiplication assembly-free. The implementation of SaRCG on multicore central processing units (CPUs) and graphics-programmable units (GPUs) is discussed, followed by numerical experiments and case-studies.