This paper presents an adaptation of the distributed-spline-based aberration reconstruction method for Shack-Hartmann (SH) slope measurements to extremely large-scale adaptive optics systems and the execution on graphics processing units (GPUs). The introduction of a hierarchical multi-level scheme for the elimination of piston offsets between the locally computed wavefront (WF) estimates solves the piston error propagation observed for a large number of partitions with the original version. To obtain a fully distributed method for WF correction, the projection of the phase estimates is locally approximated and applied in a distributed fashion, providing stable results for low and medium actuator coupling. An implementation of the method with the parallel computing platform CUDA exploits the inherently distributed nature of the algorithm. With a standard off-the-shelf GPU, the computation of the adaptive optics correction updates is accomplished in under 1 ms for the benchmark case of a 200×200 SH array.