Quantum computers leverage quantum mechanical effects to solve complex problems exponentially faster than classical computers. Their building blocks, or 'qubits', can be realized with different technologies. Silicon spin hole qubits are one of the most promising ones, thanks to their long coherence times, potentially fast manipulations, and already matured fabrication processes, as they can be encompassed within conventional CMOS transistors. Nevertheless, the performance of spin hole qubits is still far from optimal. Hence, the availability of advanced modeling platforms is key to capturing qubits' complex physics and optimizing this technology. The standard approach to simulate spin hole qubits consists of self-consistently solving the Schrödinger and Poisson equations and producing these systems' ground-state energies and charge distributions. The core operation is solving sparse eigenvalue problems for the smallest eigenpairs. For this purpose, we developed a GPU-accelerated Rayleigh–Chebyshev subspace iteration solver. Our solver relies on custom CPU/GPU kernels written in C++/CUDA and different CUDA library calls. Performance evaluations were conducted on the ALPS supercomputer and its Grace Hopper superchips. Our implementation overcomes previous time limitations achieving a speed-up of ~17x on a single GPU over the previous CPU Krylov approach, enabling high-resolution simulations of multi-qubit structures.