Abstract: To address the problem of heavy computational tasks and lengthy processing times in gravity field inversion, a parallel matrix computation method based on multi-GPU integration with CUDA was proposed. This method achieved highly parallel dense computation, significantly reducing the time required for inverting large matrices in planetary gravity field inversion by accelerating matrix multiplication and inversion operations. The computation rate was 191 times faster than that of using CPU. Moreover, it offered high computational accuracy, with inversion precision at the level of 10^–17. Applied to the computation of GRAIL lunar gravity field inversion, the proposed method, when computing matrices of truncation orders 50 and 180 respectively reduced the processing time by 94.63% and 99.51% respectively compared to CPU-based methods. Furthermore, the method successfully computed a matrix of 900th order on the High-Performance Computing Platform at Wuhan University. The method employed in this paper can effectively reduce the time needed by traditional computing models, thereby conducive to aiding in the establishment of high-order, high-precision gravity field models.