This paper presents a novel numerical method for simulation controlled-source audio-magnetotellurics (CSAMT) and radio-magnetotellurics (CSRMT) data. These methods are widely used in mineral exploration. Interpretation of the CSAMT and CSRMT data collected over an area with the complex geology requires application of effective methods of numerical modeling capable to represent the geoelectrical model of a deposit well. In this paper, we considered an approach to 3D electromagnetic (EM) modeling based on new types of preconditioned iterative solvers for finite-difference (FD) EM simulation. The first preconditioner used fast direct inversion of the layered Earth FD matrix (Green’s function preconditioner). The other combined the first with a contraction operator transformation. To illustrate the effectiveness of the developed numerical modeling methods, a 3D resistivity model of Aleksandrovka study area in Kaluga Region, Russia, was prepared based on drilling data, AMT, and a detailed CSRMT survey. We conducted parallel EM simulation of the full CSRMT survey. Our results indicated that the developed methods can be effectively used for modeling EM responses over a realistic complex geoelectrical model for a controlled source EM survey with hundreds of receiver stations. The contraction-operator preconditioner outperformed the Green’s function preconditioner by factor of 7–10, both with respect to run-time and iteration count, and even more at higher frequencies.