The Curl-Curl equation is the foundation of time-harmonic electromagnetic (EM) problems in geophysics. The efficiency of its solution is key to EM simulations, accounting for over 95% of the computation cost in geophysical inversions for magnetotelluric or controlled-source EM problems. However, most published EM inversion codes are still central processing unit (CPU)-based and cannot utilize recent computational developments on the graphic processing units (GPUs). Based on a previously proposed divergence-free algorithm developed on CPUs, this study demonstrates the current limits of the CPU-based inversion procedure. To exploit the high throughput capability of GPUs, we propose a hybrid CPU-GPU framework to solve forward and adjoint problems required for EM inversions. The large sparse linear systems arising from the staggered-grid finite difference approximation of the Curl-Curl equation are solved with a mixed-precision Krylov subspace solver implemented on a GPU. The algorithm is then tested in EM forward and adjoint calculations, with real-world three-dimensional numerical examples. Test results show promising 30× kernel-level speed-ups over the conventional CPU algorithm. This approach may further take the complex frequency domain EM inversions onto the next, practical stage on small affordable GPU platforms.