Traditional Krylov subspace methods can be slow for electromagnetic modeling, particularly at lower frequencies. One of the commonly used remedies is to apply the divergence correction iteratively during the solution process, which requires solving an additional divergence correction equation. Alternatively, we have developed an efficient regularization technique to carry out the forward modeling of the anisotropic effect for 3D controlled-source electromagnetic data. Inside this scheme, we explicitly include the gradient of a scaled divergence correction term into the original curl-curl equation for general anisotropic conductivity structure. This inclusion leads to a significantly better-conditioned linear system of equations without changing the solution of the original system and avoids the solution of an additional equation. The correctness of the developed algorithm is examined based on a vertical transverse isotropic layered model with the semianalytical solution available. Then, we use a tilted transverse isotropic ocean canonical reservoir model and a general anisotropic 3D model to investigate the stability and efficiency of the algorithm. Numerical experiments demonstrate that the developed technique is 0.49–2.9 times faster than the standard algorithm at the considered frequencies.