# How to speed up diagonalization of a Hamiltonian matrix?

I would like to diagonalize the following Hamiltonian:

\mathcal{H} = \begin{pmatrix} H & \bar{\Delta}\\ \bar{\Delta}^\dagger & -H^*\\ \end{pmatrix}
where
H = \begin{pmatrix} H_s & H_b\\ H_b^T & H_n\\ \end{pmatrix}

using the code at the end of this post, I get a matrix \mathcal{H} that is fairly sparse:

800×800 SparseMatrixCSC{ComplexF64, Int64} with 880 stored entries:
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠠⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉
⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠙⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠰⢆⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄


I want to diagonalize this matrix for a set of \phi in H_n (see code below) and I was wondering if there were some tricks/specialized types that would help speeding up this step.

Here is the code:

δ(x,y) = ===(x,y) #kronecker delta

function Δ̄(Δ, Nn, Ns, Ny)
Nx=Nn+Ns
H = zeros(Nx*Ny, Nx*Ny)
for i in 1:Ny*Ns
H[i,i] = Δ
end
return H
end

function ℍₛ(Ns, Ny, ty, tx) #\bbH
H = zeros(Ns*Ny, Ns*Ny)
for Lx in 1:Ns*Ny
for Mi in 1:Ns
for Mj in 1:Ny
Ni = trunc((Lx-1)/Ny)+1
Nj = Lx-trunc((Lx-1)/Ny)*Ny
kh=(Mi-1)*Ny+Mj

H[Lx,kh] = -ty*δ(Ni,Mi)*(δ(Nj,Mj+1)+δ(Nj,Mj-1)) - tx*δ(Nj,Mj)*(δ(Ni,Mi+1)+δ(Ni,Mi-1))
end
end
end
return H
end

function ℍb(Ns, Nn, Ny, tx)
H = zeros(Ns*Ny,Nn*Ny)
for i in 1:Ny
H[(Ns-1)*Ny+i, i]=-tx
H[i, (Nn-1)*Ny+i]=-tx
end
return H
end

function ℍₙ(ϕ, Nn, Ny, ty, tx, t1 = 6)
H = zeros(ComplexF64, Nn*Ny, Nn*Ny);
dφ = exp(-1im*ϕ/(Nn-1)*pi);
dφ̄ = conj(dφ)
for Lx in 1:Nn*Ny
for Mi in 1:Nn
for Mj in 1:Ny
Ni=trunc((Lx-1)/Ny)+1
Nj=Lx-trunc((Lx-1)/Ny)*Ny
kh=(Mi-1)*Ny+Mj

H[Lx,kh] = - ty*δ(Ni,Mi)*(δ(Nj,Mj+1)+δ(Nj,Mj-1)) - tx*dφ*δ(Nj,Mj)*δ(Ni,Mi+1) - tx*dφ̄*δ(Nj,Mj)*δ(Ni,Mi-1) + t1*δ(Lx,kh)

end
end
end
return H
end

function ℍ(ϕ, Nn, Ny, Ns, Δ, ty, tx, t1)

Hb = ℍb(Ns, Nn, Ny, tx)

H1=[ℍₛ(Ns, Ny, ty, tx) Hb]
H2=[transpose(Hb) ℍₙ(ϕ, Nn, Ny, ty, tx, t1)]
H_hop=[H1; H2]

H3=[H_hop Δ̄(Δ, Nn, Ns, Ny)]
Hconj=[conj(Δ̄(Δ, Nn, Ns, Ny)) -conj(H_hop)]
H=[H3; Hconj]

end

#To get the output posted on the forum
sparse(ℍ(pi, 20, 10, 20, 1, 4, 4, 3))

#Diagonalization, uncomment to run
#F0 = eigen!(ℍ(pi, 20, 10, 20, 1, 4, 4, 3))


Any advice is appreciated. Thanks!