# How to speed up diagonalization of a Hamiltonian matrix?

**URL:** <https://discourse.julialang.org/t/how-to-speed-up-diagonalization-of-a-hamiltonian-matrix/80777>\
**Category:** Performance\
**Tags:** linearalgebra\
**Created:** [May 9, 2022, 7:51pm UTC](https://discourse.julialang.org/t/how-to-speed-up-diagonalization-of-a-hamiltonian-matrix/80777 "2022-05-09T19:51:50Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![BastienD](https://avatars.discourse-cdn.com/v4/letter/b/c5a1d2/32.png) [@BastienD](https://discourse.julialang.org/u/BastienD)\
**Post date:** [May 9, 2022, 7:51pm UTC](https://discourse.julialang.org/t/how-to-speed-up-diagonalization-of-a-hamiltonian-matrix/80777/1 "2022-05-09T19:51:50Z")

</div>

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:

```julia
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:

```julia
δ(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!