  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 $2m$th ($m = 1, 2$) 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 $c(H)(1-C\frac{\delta^{2m-1}}{H^{2m-1}})^{2}$, where $H$ is the diameter of
subdomains and $\delta$ is the overlapping size among subdomains. The constant
$C$ is independent of the mesh size $h$ and the internal gaps among the target
eigenvalues, demonstrating that our method is optimal and cluster robust.
Meanwhile, the $H$-dependent constant $c(H)$ decreases monotonically to $1$, as
$H \to 0$, which means that more subdomains lead to the better convergence
rate. Numerical results supporting our theory are given.
