  In this work, we develop algebraic solvers for linear systems arising from
the discretization of second-order elliptic problems by saddle-point mixed
finite element methods of arbitrary polynomial degree $p \ge 0$. We present a
multigrid and a two-level domain decomposition approach in two or three space
dimensions, which are steered by their respective a posteriori estimators of
the algebraic error. First, we extend the results of [A. Mira\c{c}i, J.
Pape\v{z}, and M. Vohral\'ik, SIAM J. Sci. Comput. 43 (2021), S117--S145] to
the mixed finite element setting. Extending the multigrid procedure itself is
rather natural. To obtain analogous theoretical results, however, a multilevel
stable decomposition of the velocity space is needed. In two space dimensions,
we can treat the velocity space as the curl of a stream-function space, for
which the previous results apply. In three space dimensions, we design a novel
stable decomposition by combining a one-level high-order local stable
decomposition of [Chaumont-Frelet and Vohral\'ik, SIAM J. Numer. Anal. 61
(2023), 1783--1818] and a multilevel lowest-order stable decomposition of
[Hiptmair, Wu, and Zheng, Numer. Math. Theory Methods Appl. 5 (2012),
297--332]. This allows us to prove that our multigrid solver contracts the
algebraic error at each iteration and, simultaneously, that the associated a
posteriori estimator is efficient. A $p$-robust contraction is shown in two
space dimensions. Next, we use this multilevel methodology to define a
two-level domain decomposition method where the subdomains consist of
overlapping patches of coarse-level elements sharing a common coarse-level
vertex. We again establish a ($p$-robust) contraction of the solver and
efficiency of the a posteriori estimator. Numerical results presented both for
the multigrid approach and the domain decomposition method confirm the
theoretical findings.
