@btime changes the result

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:
Without @btime

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

        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)
        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
            Energy_bond[i] = Energy[i]/2

        #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


#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])

@btime may (will) evaluate the call multiple times, but you are setting Random.seed!(1) only once.

you may want to move the rng seeding into the setup block of @btime


I didn’t expect the random vector used in the lanczos algorithm as initial guess makes that much of a difference. Thanks