## 1 Introduction  
The 3D inversion of gravity data can determine the geometric location of subterranean anomalous bodies, and can quantitatively calculate their shape, volume, size and density distribution values, providing a strong basis for subsequent geological interpretation [1; 2; 3; 4]. The inversion of gravity data is often ill-conditioned, and Tikhonov regularization theory is typically used to reduce the multiple solutions of inversion and ensure the uniqueness and stability of inversion results. Under the regularization theory, a lot of work has been completed to address the issue of non-uniqueness in inversion. [PERSON] and [PERSON] [5] introduced the minimum volume condition on the anomalous body as a means to reduce non-uniqueness in the inversion and named it the compact inversion method. [PERSON] and [PERSON] [1] introduced the roughness matrix and a depth weighting function into the objective function of gravity inversion to effectively address the skin effect problem and named this method smooth inversion. [PERSON] and [PERSON] [6] introduced the minimum gradient support function, based on the compact inversion process, to obtain the distribution of anomalous bodies with sharp boundaries. This method was named focused inversion and was further improved in subsequent years [7; 8; 9]. [PERSON] and [PERSON] [10] examined and discussed the setting of mesh size in the inversion process by analyzing the kernel function and applied constraints such as minimum distance,flatness and smoothness to gravity data inversion through the use of the Lagrange formula. [PERSON] [11] organized the methods and techniques of 3D physical property inversion for gravity and magnetic data and proposed a magnetic data inversion method based on sparsity constraints. Smooth inversion is generally used to address the inversion problem of geological bodies with smooth boundaries, while focused inversion has been demonstrated to have improved results in restoring the 3D physical properties of geological bodies with sharp boundaries. Both of these two methods have been widely used in gravity and magnetic data inversion. The 3D inversion algorithm for gravity data can produce a reasonable density distribution and provide reference for further geological interpretation. Based on the development of fast mobile platform exploration technology (airborne and shipborne), the volume of gravity data is increasing and the accuracy of measurements are becoming more reliable. However, as the amount of observed data and model parameters increases, the inversion process faces problems of high memory consumption and long calculation time. Therefore, it is necessary to develop fast algorithms that can handle large-scale data. Using more effective optimization techniques [12; 13; 14] or compression methods [15; 16; 10] to extract more relevant information from large data sets is essential for improving inversion efficiency. Given the relative maturity of preconditioned decomposition technology, the preconditioned conjugate gradient method (PCG) has been effectively used and has produced favorable results in geophysical applications [17; 18; 19]. These methods can decrease the number of inversion iterations, thus reducing calculation time and enhancing inversion efficiency. Moreover, CUDA is a general parallel computing architecture introduced by NVIDIA, which has become a powerful tool for performing large-scale data parallel computing. Thus, CUDA has become a widely adopted tool for parallel computing particularly in the field of fast inversion of gravity or magnetic data [20; 21; 22; 23; 24]. Based on the traditional 3D inversion algorithm for gravity data that uses smooth constraints, a specific optimization scheme approach is proposed to address the low efficiency of the traditional conjugate gradient algorithm in processing large-scale data, and the application effect of the proposed method is comparable to that of the existing matrix compression methods.