In this paper we develop a parallel multilevel domain decomposition method for large-scale source identification problems governed by elliptic equations. A popular approach is to formulate the inverse problem as a PDE-constrained optimization problem. The minima satisfies a Karush–Kuhn–Tucker (KKT) system consisting of the state, adjoint and source equations which is rather difficult to solve on parallel computers. We propose and study a parallel method that decomposes the optimization problem on the global domain into subproblems on overlapping subdomains, each subdomain is further decomposed to form an additive Schwarz preconditioner for solving these smaller subproblems simultaneously with a preconditioned Krylov subspace method. For each subproblem, the overlapping part of the solution is discarded and the remaining non-overlapping part of the solution is put together to obtain an approximated global solution to the inverse problem. Since all the subproblems are solved independently, the multilevel domain decomposition method has the advantage of higher degree of parallelism. Numerical experiments show that the algorithm is accurate in terms of the reconstruction error and has reasonably good speedup in terms of the computing time. The efficiency and robustness of the proposed approach on a parallel computer with more than 1, 000 processors are reported.