The following julia script is a literal translation of the original matlab code,
using LinearAlgebra, SparseArrays
using Arpack, Plots
function bhm_basis(N, M, dim, weight)
basis = zeros(UInt8, M, dim)
basis[1, 1] = N
s = 1
while basis[M, s] < N
s1 = M - 1
while basis[s1, s] == 0
s1 = s1 - 1
end
basis[1:s1-1, s+1] = basis[1:s1-1, s]
basis[s1, s+1] = basis[s1, s] - 1
basis[s1+1, s+1] = N - sum(basis[1:s1, s+1])
s = s + 1
end
table = (weight * basis)'
ind = sortperm(table)
table = table[ind]
mini = minimum(table[2:end] - table[1:end-1]) / 4.0
return basis, table, ind, mini
end
function hint(dim, basis, M)
COO1 = zeros(3, dim)
for s1 = 1:dim
COO1[1, s1] = s1
COO1[2, s1] = s1
COO1[3, s1] = 0.5 * sum(basis[:, s1] .^ 2 - basis[:, s1])
end
return Interaction = sparse(COO1[1, :], COO1[2, :], COO1[3, :], dim, dim)
end
function hkin(dim, basis, M, table, ind, mini)
COO2 = zeros(3, dim * M * 2)
s = 0
target = [Array(2:M); 1]
for s1 = 1:dim
if mod(s1, 10000) == 0
println(s1)
end
for s2 = 1:M
if basis[s2, s1] > 0
s = s + 1
final = basis[:, s1]
final[s2] = final[s2] - 1
s3 = target[s2]
final[s3] = final[s3] + 1
value = weight * final
index = searchsortedfirst(table, value - mini)
COO2[1, s] = s1
COO2[2, s] = ind[index]
COO2[3, s] = -sqrt(final[s3] * (final[s2] + 1))
end
end
end
COO2 = COO2[:, 1:s]
Kin = sparse(COO2[1, :], COO2[2, :], COO2[3, :], dim, dim)
return Kin = Kin + Kin'
end
function gstate(Kin, Interaction, Vlist)
Cmax = zeros(length(Vlist))
ge = zeros(length(Vlist))
for s1 = 1:length(Vlist)
print(s1, '\n')
H = Kin + Vlist[s1] * Interaction
d, gstate = eigs(H, nev=1, which=:SR) # to solve the smallest algebraic eigenvalues
gstate2 = abs2.(gstate)
Cmax[s1] = maximum(gstate2)
ge[s1] = d[1]
end
return ge, Cmax
end
@time begin
N = 12
M = 12
dim = binomial(N + M - 1, M - 1)
weight = sqrt.(100 * collect(1:M) .+ 3.0)'
basis, table, ind, mini = bhm_basis(N, M, dim, weight)
Interaction = hint(dim, basis, M)
Kin = hkin(dim, basis, M, table, ind, mini)
Vlist = 0:0.5:20
ge, Cmax = gstate(Kin, Interaction, Vlist)
end
plot(Vlist, Cmax)
plot(Vlist, ge)
It is slower than matlab. Actually, it is much faster than matlab for all functions except for the function ‘gstate’, which calls the ‘eigs’ function from Arpack.
It seems that matlab is really good at matrix computation. I checked it. On my desktop, with matlab, each eigs call costs 9 seconds, while for julia, it is 14-18 seconds.
