How to speed up diagonalization of a Hamiltonian matrix?

General Usage Performance
1

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!