Effectivelly using Compile for calculate a Unitary transformation

I am new to Mathematica, and this is my first post, so if my question is not clear enough, I would be glad to read the comments and edit my question to add more information.

The problem

I need to calculate the unitary matrix $U$ that diagonalizes a SparseArray $iA$ and it columns are sorted by the eigenvalues such that [ by $^dagger$ I mean ConjugateTranspose ]:

$$ U iA U^dagger = left(
begin{array}{ccccccc}
E_1 & 0 & 0 & 0 & 0& 0
0 & E_2 & 0 & 0 & 0 & 0
0 & 0 & ddots & 0 & 0& 0
0 & 0 & 0 & -E_1 & 0 & 0
0 & 0 & 0 & 0 & -E_2 & 0
0 & 0 & 0 & 0& 0 & ddots end{array} right) $$

I made a code that worked fine and spent about 0.5 seconds for a matrix of size $1004times1004$. In the code I define my matrix $iA(k,n,J,kappa,h)$ of size $(2n+4)times(2n+4)$, then I use Eigensystem and apply mySort to sort the eigenvectors in the way I wanted and formed my unitary matrix with this eigenvectors. I need it to be faster because I want to examine how the coefficients of the unitary matrix $U$ depend upon the parameters in the interval ${ 0leq k < pi , , , 0< kappa leq 1 , , , 0< hleq 1 }$, for $J=1$ and optimally for large $n$. In practice, I would like to test for about $100$ values of $k$, $10$ for $kappa$ and $10$ values for $h$; so I will need to calculate about $10.000$ of matrices.

• Q1: There is a more efficient way to calculate this Unitary matrix?

The first problem is that the Mathematica automatically convert my SparseArray, to a normal one, since I want to calculate all the eigenvectors (which is needed). Return the standard warning:

Eigensystem::arh: Because finding 1004 out of the 1004 eigenvalues and/or eigenvectors is likely to be faster with dense matrix methods, the sparse input matrix will be converted...

I try to use other Methods for the Eigensystem, but neither of this is optimized for calculating all eigenvector. From what I have seen in others post this would be the maximum performance if I indeed wanted all eigenvectors, then I try to reduce the numbers of eigenvectors that I want. Given the symmetry of the matrix $iA$ I also tried to calculate just half of the eigenvectors, corresponding to only positive eigenvalues, using Method->{"FEAST","Interval"->{0.,10.}} but I do not see any time reduction.

• Q2: There is a more efficient method for calculating all eigenvectors corresponding to positive eigenvalues?

I also want to put this whole code inside a function that takes values $(k,n,J,kappa,h)$ and return the unitary matrix I have tried to use the Compile, but the Compile return some errors related with the Type of the matrix

Compile::cptype: N is not supported for type None; evaluation will use the uncompiled function.

and

Compile::cplist: SortBy[Compile`FunctionVariable$7579,First] should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function.

I have read the documentation for Compile, and I don’t know if it is appropriate for this kind of problem, but I was thinking that it would speed up the calculation. From what I have seen, the Compile can not return a matrix as an output (it is true?), but for me would be enough to analyse element by element of the matrix.

• Q3: Can I use Compile in this example? Or at least use it to generate a compiled function that outputs the $ij$-th element ?

My code

As I described above, my code has a few fundamental elements: the matrix itself, the sorting procedure and the code to calculate the unitary matrix.

• The Matrix

    iA[k_, n_, J_, κ_, h_] :=   SparseArray[{

{i_, i_} ->      If[1 < i < 2*n + 4,  If[OddQ[i], -2*κ*Sin[k], 2*κ*Sin[k]], 0], 
{i_, j_} /; Abs[i - j] == 1 ->  If[(i == 1 || j == 1) || (i == 2*n + 4 || j == 2*n + 4), If[i == 1 || i == 2*n + 4, (-I)*h, I*h],  If[EvenQ[Min[i, j]],  2*J*Cos[k/2]*(i - j)*I, (-(i - j))*I*J]], 
{i_, j_} /;  Abs[i - j] == 2 -> If[(i == 1 || j == 1) || (i == 2*n + 4 || j == 2*n + 4), 0, If[OddQ[i], 2*κ*Sin[k/2], -2*κ*Sin[k/2]]]}

, {2*n + 4, 2*n + 4}]; 

• MySort

mySort = (Transpose[ Chop[Flatten[{Drop[SortBy[#1, First], Length[#1]/2],Reverse[Take[SortBy[#1, First], Length[#1]/2]]}, 1]]] & )[Transpose[#1]] & ; 

• The code for the calculate the Unitary matrix
  The code that works but is slow:

U[k_, n_, J_, κ_, h_] = mySort[Eigensystem[N[iA[k, n, J, κ, h][[2]]]]]

and the code with Compile that don’t even work

U = Compile[  { {k, _Real, 0}, {n, _Integer, 0}, {J, _Real, 0}, {κ, _Real, 0}, {h, _Real, 0} , {m, _Integer, 0} , {y, _Integer, 0} },  mySort@Eigensystem@N@Normal@ iA[k, n, J, κ, h]  , CompilationOptions -> {"InlineExternalDefinitions" -> True}];

Additional Info

• The matrix $iA$ is a Hermitian and pentadiagonal of the form [ it has only three distinct $2times 2$ blocks the repeates across the matrix] :

! $ iA = left(
begin{array}{cccccccc}
0 & -i h & 0 & 0 & 0 & 0 & 0 & 0
i h & 2 kappa sin (k) & -2 i J cos left(frac{k}{2}right) & -2 kappa sin left(frac{k}{2}right) & 0 & 0 & 0 & 0
0 & 2 i J cos left(frac{k}{2}right) & -2 kappa sin (k) & i J & 2 kappa sin left(frac{k}{2}right) & 0 & 0 & 0
0 & -2 kappa sin left(frac{k}{2}right) & -i J & 2 kappa sin (k) & -2 i J cos left(frac{k}{2}right) & -2 kappa sin left(frac{k}{2}right) & 0 & 0
0 & 0 & 2 kappa sin left(frac{k}{2}right) & 2 i J cos left(frac{k}{2}right) & -2 kappa sin (k) & i J & 2 kappa sin left(frac{k}{2}right) & 0
0 & 0 & 0 & -2 kappa sin left(frac{k}{2}right) & -i J & 2 kappa sin (k) & -2 i J cos left(frac{k}{2}right) & 0
0 & 0 & 0 & 0 & 2 kappa sin left(frac{k}{2}right) & 2 i J cos left(frac{k}{2}right) & -2 kappa sin (k) & i h
0 & 0 & 0 & 0 & 0 & 0 & -i h & 0
end{array}
right) $

• I also don’t know if the way I defined the $iA$ matrix and $mySort$ is the most efficient, but at least I have tested I lot and they work.
• The spectrum of $iA$ has a energy gap for all $k$ but for $k=0$.
• The matrix $iA$ has a symmetric spectrum of eigenvalues (i.e. if $E$ is a eigenvalue, then $-E$ will also be).
• I do not know a relation between the eigenvector corresponding to the eigenvalues $E$ e $-E$, but if someone has any hint I would appreciate a lot!
• If it is not clear, I want to calculate the coefficients of $U$ as a function of $(k,kappa,h)$ for fixed $n$ and $J$. I want to see which coefficients are more relevant [biggest absolute value] in the regime of $k to 0$. Then I will plot how those $U_{i,j}$ depends upon each parameter, and hopefully I with figure an analytic dependence in that regime ( expect a linear behaviour in $k$, but no hint in the $kappa$ and $h$ dependence).