We consider the problem of accelerating the iterative diagonalization of Hamiltonian operators for electronic structure calculations in plane-wave Density Functional Theory. The complexity bottleneck of existing subspace iteration schemes is that the Rayleigh-Ritz procedure for extracting eigenvectors from a subspace has limited scalability with the number of processors. On the other hand, polynomial filtering methods for constructing subspaces can greatly benefit from the parallel efficiency of Hamiltonian applications using batch processing. The latter feature becomes even more advantageous on Graphics Processing Unit architectures in particular. For this reason, we focus on the spectrum slicing method, which is a special type of polynomial filtering allowing to apply the Rayleigh-Ritz step to fewer vectors. We present a new implementation of this scheme in the ABINIT plane-wave code using GPUs. We also propose a complexity model that allows to achieve load balance in the number of filtering operations and the number of vectors on which we apply the Rayleigh-Ritz step per slice. The numerical performance of our implementation is compared to that of existing algorithms in ABINIT, in particular Chebyshev Filtering and Locally Optimal Block Preconditioned Conjugate Gradient, for systems of up to ten thousand bands.