We present a fast inversion algorithm tailored for high-resolution potential field (gravity or magnetic) anomaly maps. The algorithm design objectives are to optimize performance with respect to the size of problem that can be solved and computation speed. It is based on a finite element method (FEM) discretization of both the inversion and forward problems using unstructured tetrahedral meshes with coarsening at the far field. This allows a significant reduction in computational costs from a structured mesh with the same ground level resolution, as well as supporting geometrical features such as topography. Minimization of the cost function uses the Lagrange minimization approach and the resulting system of partial differential equations is solved using the integral preconditioned conjugate gradient method (I-PCG). I-PCG works naturally with an unstructured mesh FEM and can effectively precondition with the Hessian of the regularization term. In contrast to established approaches, minimization of the cost function occurs prior to discretization and using I-PCG avoids the need to invert a large dense sensitivity matrix, a major obstacle to solving large 3D problems using parallel computers. Fast solvers for the forward problem, the adjoint problem and the preconditioner make use of the smoothed aggregation algebraic multi-grid method for the preconditioner in another PCG iteration (AMG-PCG). Parallel implementation uses a domain decomposition approach to support scalability with spatial extent. The proposed method is applied to a synthetic example and tohigh-resolution data set of rastered gravity and magnetic anomaly maps from the Eastern Goldfields in Western Australia with the finest resolution of over 1 million data points and about 60 million cells for the 3D inversion. For this field application we demonstrate that the proposed algorithm achieves computation time scaling with grid size, an indication of optimal algorithm performance and delivers strong scalability with number of cores. The gravity inversion result is compared to the inversion result using the Broyden Fletcher Goldfarb Shanno (BFGS) method on a structured grid.