  In this paper, we develop an EXCMG method to solve the three-dimensional
Poisson equation on rectangular domains by using the compact finite difference
(FD) method with unequal meshsizes in different coordinate directions. The
resulting linear system from compact FD discretization is solved by the
conjugate gradient (CG) method with a relative residual stopping criterion. By
combining the Richardson extrapolation and tri-quartic Lagrange interpolation
for the numerical solutions from two-level of grids (current and previous
grids), we are able to produce an extremely accurate approximation of the
actual numerical solution on the next finer grid, which can greatly reduce the
number of relaxation sweeps needed. Additionally, a simple method based on the
midpoint extrapolation formula is used for the fourth-order FD solutions on
two-level of grids to achieve sixth-order accuracy on the entire fine grid
cheaply and directly. The gradient of the numerical solution can also be easily
obtained through solving a series of tridiagonal linear systems resulting from
the fourth-order compact FD discretizations. Numerical results show that our
EXCMG method is much more efficient than the classical V-cycle and W-cycle
multigrid methods. Moreover, only few CG iterations are required on the finest
grid to achieve full fourth-order accuracy in both the $L^2$-norm and
$L^{\infty}$-norm for the solution and its gradient when the exact solution
belongs to $C^6$. Finally, numerical result shows that our EXCMG method is
still effective when the exact solution has a lower regularity, which widens
the scope of applicability of our EXCMG method.
