We describe solution algorithms based on block preconditioning for discretizations of the fully coupled steady-state magnetohydrodynamics (MHD) equations. The equations are posed in an exact-penalty formulation, which has the advantage of implicitly enforcing the divergence-free condition of the magnetic field without a Lagrange multiplier. The solvers are constructed by generalizing block methods developed for the Navier-Stokes equations to handle models of MHD, using ideas of "approximate commutators" in conjunction with automated strategies designed to compensate for errors in commutators. We show that the resulting preconditioned Krylov subspace solvers are robust for a wide range of physical parameters and with respect to mesh refinement.