The use of an incomplete iterative solver aimed at approximating solution of equilibrium equation resulting from the "nite element semidiscretization in space is studied. Our approach is motivated by the fact that: (i) the e!ective sti!ness matrix is well conditioned because of the stabilizing e!ect of the diagonally dominated mass matrix, (ii) the dominance of the temporal error in implicit computations with large step size, and (iii) the utilization of a class of iterative methods where the primary computational cost is associated with a construction of the preconditioner rather than with an iterative process. One of the primary goals of the present manuscript is to control the temporal and the equilibrium solution errors and to give their a posteriori estimates. Numerical experiments reveal that the computational cost of three cycles conducted with the Generalized Aggregation Multilevel solver [1}4] is comparable to a single back substitution on the source grid, while the resulting equilibrium solution error is negligible compared with the local temporal error.