A three-dimensional (3D) inversion algorithm for magnetotelluric (MT) data has been developed based on the Gauss-Newton method. The algorithm employs a 3D MT modeling algorithm based on an edge finite element method (FEM) with hexahedral elements being able to delineate surface topography. The modeling algorithm solves the resulting system of equations using a bi-conjugate gradient (BICG) method with a symmetric Jacobian preconditioner. To verify the inversion algorithm, we make inversion of MT data for a flat-earth surface model with a 3D conductive body in a homogeneous host. Further numerical tests are made for MT data over a surface topography. The model under consideration has a conductive or resistive dike under a 3D trapezoidal hill. For inversion of distorted MT data by surface topography, we first remove the topography effects using an impedance correction method. Numerical results indicate that inversion of MT data over the flat-surface model located the conductive target. Inversion of topography-effects-corrected MT data explains the presence of the conductive or resistive target while inversion of distorted MT data gives no clue.