In this thesis we present a multigrid approach for systems of linear equations involving the Wilson Dirac operator from lattice quantum chromodynamics (lattice QCD), the fundamental physical theory of the strong interaction between the quarks. The arising systems are non-Hermitian, large and ill-conditioned. Our multigrid approach combines two components that have already been used separately in lattice QCD, namely the domain decomposition method "Schwarz Alternating Procedure" as a smoother and the aggregation based interpolation. This approach outperforms the most popular Wilson solvers for lattice QCD and stimulated the implementation of AMG solvers in other QCD computation tools. The thesis is structured as follows: First, we give an introduction to the construction of the Wilson Dirac operator and point out its mathematical properties. We discuss various kinds of domain decomposition methods and give implementation details for a parallel implementation, followed by numerical results from our MPI-C implementation. Then we give an introduction to multigrid methods for lattice QCD and in particular to our domain decomposition type approach including its adaptive setup procedure and implementation details. The differences to other existing hierarchical approaches for lattice QCD are pointed out in detail, theoretically as well as numerically with tests from MPI-C implementations. Additionally, we present several scaling studies for various parameters and large lattice configurations on a parallel cluster computer with up to 8,192 cores. As a non-standard preconditioning technique we apply our Wilson solver as an auxiliary space type preconditioner to the overlap operator from lattice QCD. As opposed to the Wilson Dirac operator, the overlap operator preserves the important physical property of chiral symmetry at the expense of requiring much more effort when solving systems with this operator. We give a mathematical analysis which shows that our preconditioner is effective in an idealized setting where the Wilson Dirac operator is assumed to be normal. This is then confirmed by numerical experiments which are performed for large lattice configurations coming from state-of-the-art physical simulations. Compared to standard Krylov subspace methods we obtain a speed up in runtime of more than an order of magnitude.