## 1 Introduction  
TEM has undergone rapid development and achieved significant progress, becoming widely used in various fields such as mineral exploration [1, 2, 3, 4], hydrogeophysical exploration [5, 6, 7], tunnel prediction [8, 9], and more [10, 11]. TEM forward modeling serves as the foundation for data interpretation and 3D inversion. As the demand for detailed detection and data interpretation continues to grow, the need for 3D large-scale forward modeling has become increasingly urgent [12]. Currently, TEM 3D forward modeling methods can be broadly categorized into three types: the finite difference method [13, 14, 15, 16], the implicit time-step method [17, 18, 19, 20, 21], and the Krylov subspace method [22, 23, 24, 25, 26, 27]. Among these methods, the implicit time-step and Krylov subspace methods typically employ direct solvers to simulate the TEM response for small- to medium-scale models. However, simulating the electromagnetic response of large-scale models presents considerable challenges in terms of memory usage [28]. While the iterative method based on the implicit time-step can also address large-scale models, its computational efficiency requires further improvement, and selecting the optimal time step remains a challenging task [17; 19; 29]. The three-dimensional finite difference method in the time domain does not require solving linear equations, but careful selection of the time step and mesh size is still needed to prevent numerical overflow due to numerical dispersion [14; 15; 16]. In recent years, the Krylov subspace method has made significant advancements compared to the traditional implicit time-step method. [PERSON] and [PERSON] addressed the slow convergence issue of the SAI-Krylov subspace algorithm, proving that the Krylov subspace method always approximates an exponential function with a superlinear convergence rate, independent of mesh size [30; 31]. Additionally, [PERSON] et al. highlighted that the Krylov subspace method outperforms other methods in terms of convergence [25]. Furthermore, the Krylov subspace method eliminates the need to select a time step, enabling the calculation of the electromagnetic response at any point in the off-time period [25; 26; 27]. [PERSON] et al. demonstrated fast TEM forward modeling using the SAI-Krylov subspace method [26]. Compared with the rational Krylov subspace algorithm, the SAI-Krylov subspace method offers faster computational efficiency [32]. Based on these advancements, this paper proposes a new method for solving TEM responses using the SAI-Krylov subspace method. However, for large-scale models with millions to even billions of unknowns, the combined computational strategy of SAI-Krylov and direct solvers is relatively inefficient. Compared to the direct solvers, the iterative method only requires storing the sparse matrix and a small number of vector elements [17]. The primary computational effort in the iterative process involves multiplying the coefficient matrix by the vector. The PCG method is a highly effective optimization technique that only requires first-order derivative information. This method not only addresses the slow convergence of the steepest descent method but also avoids the limitations of the Newton method, which requires storing and calculating the Hessian matrix [33]. The PCG method is well suited for solving large linear equations and is also applicable to large nonlinear problems. Compared to other iterative methods, the PCG method offers clear advantages in accuracy and efficiency when using the same initial value and error threshold [34]. As a classic iterative method, the PCG technique has been widely employed in electromagnetic numerical simulations. [PERSON] and [PERSON] applied the conjugate gradient method to three-dimensional resistivity inversion, achieving satisfactory results [35].