I will write a difficult to read program because I know its analytical solution.
What you just need to to is putting @btime in front of the function call, at the end.
Read the last result of the 4th column:
With @btime:
-0.269957827874900
Without @btime
-0.269946270346980
It’s quite funny the fact that the result with @btime is more accurate than the one without it. The algorithm is a really simple and approximate variational algorithm.
using Printf
using SparseArrays
using Random
using LinearAlgebra
using BenchmarkTools
using KrylovKit
function IDMRG(Delta::Real, max_states::Integer, NIter::Integer, Energy::Vector, Energy_2::Vector, Energy_bond::Vector, lattice_size::Vector, Truncation_Error::Vector)
Sz = sparse(Float64[[0.5, 0] [0, -0.5]])
Sp = sparse(Float64[[0, 0] [1.0, 0]])
Sm = sparse(Float64[[0, 1.0] [0, 0]])
Ham = sparse(zeros(2,2))
Block_Sp = Sp
Block_Sm = Sm
Block_Sz = Sz
Block_Ham = Ham
m::Int64 = 2 #hilbert size
temp_matr = Array{Float64}(undef, max_states*2, max_states)
for i = 1:NIter
lattice_size[i] = 2*i +2 #spatial row_size
#Adding 1 site
H_new = kron(Block_Ham, sparse(I,2,2)) .- Delta*kron(Block_Sz, Sz) .+ 0.5*(kron(Block_Sp, Sm) .+ kron(Block_Sm, Sp))
Block_Sz = kron(sparse(I,m,m), Sz) #definition of new Block_Sz
Block_Sp = kron(sparse(I,m,m), Sp) #definition of new Block_Sp
Block_Sm = kron(sparse(I,m,m), Sm) #definition of new Block_Sm
#SuperBlock
H_super = kron(H_new, sparse(I,m*2,m*2)) .+ kron(sparse(I,m*2,m*2), H_new) .- Delta*kron(Block_Sz, Block_Sz) .+ 0.5*(kron(Block_Sp, Block_Sm) .+ kron(Block_Sm, Block_Sp))
droptol!(H_super, 5e-16)
#trunc index
m_next::Int64 = m*m*4 # hilbert space row_size SuperBlock
rho_dim::Int64 = m*2 #sqrt(m_next)
trunc_index::Int64 = min(rho_dim,max_states) #this is applied to density matrix diag
#Diagonalizing SuperBlock Lanczos
Energy_ground, evec_computed, info = eigsolve(H_super, rand(Float64, m_next), 1, :SR, Lanczos(krylovdim=10, tol=1e-16)) # maxiter = 100
if info.converged < 1
println("Lanczos did not converge in iteration: ", i)
exit()
end
Energy[i] = Energy_ground[1] #output of eigensolve is a vector of eigenvals
psi_ground = evec_computed[1]
Energy_2[i] = Energy[i]/lattice_size[i]
if i> 1
Energy_bond[i] = (Energy[i] - Energy[i-1]) / 2
else
Energy_bond[i] = Energy[i]/2
end
#Reduced Density Matrix
psi_matr = reshape(psi_ground, (rho_dim, rho_dim)) #psi matr is a view of psi allocated data
rho_matr = psi_matr * psi_matr' #a way of tracing rho_matrix on the right part
rho_matr_sym = Symmetric(rho_matr) #it seems this is indeed symmetric
#Diagonalizing Density Matrix, Trunc matrix
evals, evecs = eigen!(rho_matr_sym, sortby=-) #if I use eigen! i need to not use rho_matrx again
trunc_evals = @view evals[1:trunc_index]
trunc_matrx = @view evecs[:,1:trunc_index]
Truncation_Error[i] = 1 - sum(trunc_evals)
# Preallocation for mul!
temp = @view temp_matr[1:rho_dim,1:trunc_index]
#Higher dimension kron
mul!(temp, H_new, trunc_matrx)
Block_Ham = trunc_matrx' * temp
mul!(temp, Block_Sz, trunc_matrx)
Block_Sz = trunc_matrx' * temp
mul!(temp, Block_Sp, trunc_matrx)
Block_Sp = trunc_matrx' * temp
mul!(temp, Block_Sm, trunc_matrx)
Block_Sm = trunc_matrx' * temp
m = trunc_index # hilbert space row_size before next kron products
end
end
Random.seed!(1)
#Definition of starting parameters
Delta::Float64 = 0.5
max_states::Int16 = 5
NIter::Int32 = 30 #Final lattice size: 2*NIter + 2
#output vectors
lattice_size = Vector{Int64}(undef, NIter)
Energy = Vector{Float64}(undef, NIter)
Energy_2 = Vector{Float64}(undef, NIter)
Energy_bond = Vector{Float64}(undef, NIter)
Truncation_Error = Vector{Float64}(undef, NIter)
@btime IDMRG(Delta, max_states, NIter, Energy, Energy_2, Energy_bond, lattice_size, Truncation_Error)
for i = 1:NIter
print(lpad(lattice_size[i], 2, "0"), " ")
@printf("| %.15f | %.15f | %.15f | %.15f | \n", Energy[i], Energy_bond[i], Energy_2[i], Truncation_Error[i])
end