SUMMARY Geomagnetic depth sounding (GDS) is a geophysical electromagnetic (EM) method that studies the deep structure and composition of the Earth by using long-period EM signals from geomagnetic observatories and satellites. In this paper, a 3-D anisotropic GDS modelling algorithm is developed. The curl–curl equation is discretized using the edge-based finite-element method on unstructured tetrahedral grids. In order to solve the computationally demanding problem of EM modelling on a global scale, the complex linear system is first separated into the equivalent real linear systems and then the real system is iteratively solved by the flexible generalized minimum residual method with a block diagonal pre-conditioner. This will greatly reduce the condition number of the linear system and thus speed up the solution process. We verify the accuracy of the proposed algorithm by comparing our results with the existing methods. After that, we design a subduction zone model to simulate the EM responses under isotropic and anisotropic environments, respectively. The numerical results show the high efficiency of the proposed algorithm and the response differences between isotropic and anisotropic models. This research can provide theoretical and technical support for the high-accuracy and efficient inversion of GDS data for the geo-dynamic study.