We have developed a new algorithm for the inversion of magnetotelluric (MT) data. The developed algorithm is built to be fast, versatile, and accurate. A fast inversion algorithm has to include a fast forward-modeling routine. To achieve that, a hybrid approach consisting of finite-difference (FD) and finite-element (FE) methods is used to benefit from the speed of the FD method and the flexibility to add topographic features of the FE method. To reduce the number of cells, and thus reducing the size of the system to be solved in the forward and pseudoforward solutions, different meshes for various groups of frequencies are used. Then, these are mapped onto the inversion mesh by a mesh-decoupling technique to further accelerate the inversion. To build a versatile inversion algorithm, the capability of using different data types is implemented. In addition to the impedance tensor and the magnetic transfer function, the algorithm also computes the phase tensor and phase vector, which are distortion-free forms of MT data. It is also possible to invert intersite data and their respective phase tensors using the developed code. Furthermore, the distortion matrix can also be estimated as a parameter. The new code is tested with different noisy and distorted synthetic data measured on a surface with topography to evaluate the inversion accuracy and computational efficiency. The results indicate that the code is accurate and that the runtimes are reasonable for the large 3D models considered. Using four graphics processing units, the hybrid forward-modeling approach and the mesh-decoupling technique together result in a 12 times speedup for the examples presented in this study.