We apply a two-dimensional Cartesian finite element treatment to investigate infinite Prandtl number thermal convection with temperature, strain rate and yield stress dependent rheology using parameters in the range estimated for the mantles of the terrestrial planets. To handle the strong viscosity variations that arise from such nonlinear rheology in solving the momentum equation, we exploit a multigrid method based on matrix-dependent intergrid transfer and the Galerkin coarse grid approximation. We observe that the matrix-dependent transfer algorithm provides an exceptionally robust and efficient means for solving convection problems with extreme viscosity gradients. Our algorithm displays a convergence rate per multigrid cycle about five times better than what other published methods (e.g., CITCOM of Moresi and Solomatov, 1995) offer for cases with similar extreme viscosity variation. The algorithm is explained in detail in this paper.