<p>I am not entirely sure how <em>uniform sampling of the eigenvalues</em> would work, but I think you are looking for <a href="https://www.caam.rice.edu/software/ARPACK/" rel="nofollow noreferrer">ARPACK</a>. ARPACK would use matrix-vector products to find your eigenvalues, so I am not entirely sure if the Real/Im decomposition is required in this case (hard to say without knowing a lot about the <code>U</code>).</p>

<p>Also, you might want to look at <a href="http://www.ecs.umass.edu/~polizzi/feast/" rel="nofollow noreferrer">FEAST</a> algorithm, which would benefit a lot from the <a href="https://arxiv.org/abs/1203.4031" rel="nofollow noreferrer">given search contour</a>.</p>

<p>I am not aware of the existing linking of Julia to those libraries, but I don't think it is a problem since Julia can call C functions.</p>

<p>Here, I gave some brief ideas, and <a href="https://scicomp.stackexchange.com/">Computational Science</a> might be a better place to find the right crowd. However, a lot more details about <code>U</code>, its sparsity, size, and what does "uniform sampling of eigenvalues in the interval" means would be required.</p>
Homogeneous eigenvalue sampling of a sparse unitary matrix
<p>I work with Julia, but I think the question is more general. Suppose that one wants to find the spectrum of a very large (sparse) unitary matrix <code>U</code> numerically. As is reported in many entries, diagonalizing by brute force using <code>eigs</code> ends without eigenvalue convergence.</p>

<p>The trick would be then to work with simpler expressions, i.e. with </p>

<pre><code>U_Re = real(U + U')*0.5                                                                             
U_Im = real((U - U')*-0.5im)
</code></pre>

<p>My question is, is there a way to obtain a uniform sampling in finding the eigenvalues? That is, I would like to obtain, say <code>10e3</code> eigenvalues for <code>U_Re</code> and <code>U_Im</code> in the interval <code>[-1,1]</code>.</p>