The novel contribution of this paper relies in the proposal of a fully
implicit numerical method designed for nonlinear degenerate parabolic
equations, in its convergence/stability analysis, and in the study of the
related computational cost. In fact, due to the nonlinear nature of the
underlying mathematical model, the use of a fixed point scheme is required and
every step implies the solution of large, locally structured, linear systems. A
special effort is devoted to the spectral analysis of the relevant matrices and
to the design of appropriate iterative or multi-iterative solvers, with special
attention to preconditioned Krylov methods and to multigrid procedures: in
particular we investigate the mutual benefit of combining in various ways
suitable preconditioners with V-cycle algorithms. Numerical experiments in one
and two spatial dimensions for the validation of our multi-facet analysis
complement this contribution.