In this paper, we propose a two-level block preconditioned Jacobi–Davidson (BPJD) method for efficiently solving discrete eigenvalue problems resulting from finite element approximations of th () order symmetric elliptic eigenvalue problems. Our method works effectively to compute the first several eigenpairs, including both multiple and clustered eigenvalues with corresponding eigenfunctions, particularly. The method is highly parallelizable by constructing a new and efficient preconditioner using an overlapping domain decomposition (DD). It only requires computing a couple of small scale parallel subproblems and a quite small scale eigenvalue problem per iteration. Our theoretical analysis reveals that the convergence rate of the method is bounded by , where is the diameter of subdomains and is the overlapping size among subdomains. The constant is independent of the mesh size and the internal gaps among the target eigenvalues, demonstrating that our method is optimal and cluster robust. Meanwhile, the -dependent constant decreases monotonically to 1, as , which means that more subdomains lead to the better convergence rate. Numerical results supporting our theory are given.