Quantitative interpretation of large-scale controlled-source electromagnetic (CSEM) data in the frequency domain requires efficient and stable three-dimensional (3D) forward modeling and inversion codes. In this paper, we present a 3D forward modeling scheme for frequency-domain CSEM surveys based on the mimetic finite volume (MFV) method, which solves the Helmholtz equation for the total electric field. An important step in this 3D modeling scheme is to solve the large linear equation system resulting from MFV discretization. In order to obtain stable and accurate numerical solutions, the linear equation system is solved by a direct-matrix solver, namely MUMPS, which is more robust than commonly used iterative solvers (e.g, Krylov subspace iterative techniques) for numerically difficult cases. The algorithm is very suitable for multi-source CSEM modeling by separating the solving process into a single expensive matrix factorization and relatively inexpensive forward backward substitutions for many right-hand sides. For total-field formulation, dense gridding in the vicinity of source points is usually required to mitigate source singularities, and extra padding cells at the boundaries are needed to meet boundary conditions. Both of them can quickly increase the size of the linear equation system to be solved, making it computationally expensive for direct solutions. To avoid substantial increase of the size of the discretized model, a source correction technique is applied to reduce source singularities and boundary effects. In addition, considering the independence of computation for different frequencies, a parallel forward modeling scheme based on frequency partition is implemented to speed up the simulation of multi-frequency CSEM problems. Numerical experiments have been carried out to evaluate the performance of our forward modeling algorithm. Numerical solutions using this algorithm show good agreement with quasi analytic solutions for 1D layered models, and numerical errors are significantly reduced with source correction in the vicinity of source points. In addition, comparison of simulated data generated by our algorithm to published 3D data for a typical marine 3D model validates our algorithm. Besides, statistical results demonstrate that parallel computing based on simple frequency partition approach can achieve nearly linear speedup due to the independence of computation between frequencies for multi-frequency CSEM simulation.