Summary: Preconditioners based on various multilevel extensions of two-level finite element methods lead to iterative methods which have an optimal order computational complexity with respect to the size (or discretization parameter) of the system. The methods can be on block matrix factorized form, recursively extended via certain matrix polynomial approximations of the arising Schur complement matrices or on additive, i.e., block diagonal form using stabilizations of the condition number at certain levels. The resulting spectral equivalence holds uniformly with respect to jumps in the coefficients of the differential operator and for arbitrary triangulations. An important part of the algorithm is the treatment of the systems with the diagonal block matrix, which arises on each finer level in a recursive refinement method and corresponds to the added degrees of freedom on that level. This block is well-conditioned for model type problems but becomes increasingly ill-conditioned when the coefficient matrix becomes more anisotropic or, equivalently, when the mesh aspect ratio increases. In the paper some methods are presented to approximate this matrix also leading to a preconditioner with spectral equivalence bounds which hold uniformly with respect to both the problem and discretization parameters. The same holds therefore also for the preconditioner to the global matrix. Such uniform bounds have not been achieved by other methods.