We propose a fast and accurate 3D large-scale forward modeling method to compute the gravity gradient tensor, in terms of a fourth-order compact finite difference scheme and the extrapolation cascadic multigrid (EXCMG) method. Firstly, the 19-point fourth-order compact finite difference scheme with unequal mesh sizes in different coordinate directions is used to discretize the governing 3D Poisson equation. The resulted symmetric positive definite linear systems are solved by the EXCMG method. Then, the second-order derivatives (gravity gradient tensor) are calculated by solving a series of tridiagonal linear systems resulting from the fourth-order compact finite difference discretization. Finally, numerical result for a cubic model with positive density is used to verify the accuracy of our presented method. It shows that our method can give nearly fourth-order accurate approximations to gravitational potential and gravity gradient tensor. The accuracy is much higher and more efficient than traditional finite difference methods.