MKL package wrong results

General Usage
, , ,
1

I have this program. I post it as it is.
If I write:

using MKL
using MKLSparse

It gives me (that’s an iteration of results with a matrix parameter on the first column):

02 | -0.37499999999999944 | 0.0 | -0.18749999999999972
04 | -0.375 | -2.7755575615628914e-16 | -0.09375
08 | -0.875 | -0.125 | -0.109375
16 | -3.2757653911111717 | -0.30009567388889646 | -0.20473533694444823
32 | -7.198931123379452 | -0.24519785826676754 | -0.2249665976056079

If i don’t use them, it gives me:

02 | -0.37499999999999944 | 0.0 | -0.18749999999999972
04 | -0.9098524898764077 | -0.2674262449382041 | -0.2274631224691019
08 | -2.68529096887001 | -0.4438596197484006 | -0.33566137110875127
16 | -5.43894887951664 | -0.3442072388308287 | -0.33993430496979
32 | -11.008233268476827 | -0.3480802743100117 | -0.34400728963990085

If instead of kron, you add kronecker package, nothing changes without MKL.
If you add MKL too, then it gives:

02 | -0.37499999999999944 | 0.0 | -0.18749999999999972
04 | -0.7838137951917248 | -0.20440689759586267 | -0.1959534487979312
08 | -1.603161567518434 | -0.20483694308167733 | -0.20039519593980426
16 | -3.9978146703099906 | -0.29933163784894457 | -0.2498634168943744
32 | -8.178459637176488 | -0.26129031042915607 | -0.25557686366176524

In the following text i will provide 2 versions of the code. The first one allocates a temporary matrix each iteration, called temp. See around line 45.
The first version gives the results above.
The second version prallocates temp with the maximum size it reaches through all the iterations. I view inside of it each iteration.
Without using MKL it works and the results are the same.
When using MKL it gives:

ERROR: LoadError: InexactError: Float64(-0.9999999999999998 + 1.4994064124073605e-43im)
Stacktrace:
 [1] Real
   @ ./complex.jl:44 [inlined]
 [2] convert
   @ ./number.jl:7 [inlined]
 [3] setindex!(A::Vector{Float64}, x::ComplexF64, i1::Int64)
   @ Base ./array.jl:1021
 [4] NRG(Delta::Float64, max_states::Int16, NIter::Int32, Energy::Vector{Float64}, Energy_2::Vector{Float64}, Energy_bond::Vector{Float64})
   @ Main ~/Desktop/Programmi/MODULO_05/executable_and_libraries/bin/simple_nrg_optimized.jl:36
 [5] top-level scope

Instead of writing two different version, i will write the additional 2 lines of the second version and comment on the right

using SparseArrays
using LinearAlgebra
using BenchmarkTools
using MKL
using MKLSparse

function NRG(Delta::Real, max_states::Integer, NIter::Integer, Energy::Vector, Energy_2::Vector, Energy_bond::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

    #temp_matr = Array{Float64}(undef, max_states*max_states, max_states) # ENABLE IN VERSION 2

    m::Int64 = 2 #hilbert size

    for i = 1:NIter
        lattice_size::Int64 = 2^i #spatial row_size

        #Hamiltonian constructor
        H_super = kron(Block_Ham, sparse(I,m,m)) .+ kron(sparse(I,m,m), Block_Ham) .- Delta*kron(Block_Sz, Block_Sz) .+ 0.5*(kron(Block_Sp, Block_Sm) .+ kron(Block_Sm, Block_Sp))
        H_dense_super = collect(H_super) #I make a dense matrix as this is a dummy example with full diagonalization.

        #trunc index
        m_next = m*m # H_super size
        trunc_index = min(m_next,max_states)

        #Diagonalization part
        evals, evecs = eigen!(H_dense_super) #ascendent ordered
        Energy[i] = evals[1]
        trunc_matrx = @view evecs[:,1:trunc_index]

        Energy_2[i] = Energy[i]/lattice_size
        if i> 1
            Energy_bond[i] = (Energy[i] - Energy[i-1]) / (lattice_size/2)
        end

        # Preallocation for mul!
        #temp = @view temp_matr[1:m_next,1:trunc_index] # ENABLE IN VERSION 2
        temp = Array{Float64}(undef, m_next, trunc_index) # DISABLE IN VERSION 2

        #Higher dimension kron
        mul!(temp, H_super, trunc_matrx)
        Block_Ham = trunc_matrx' * temp

        mul!(temp, kron(sparse(I, m, m), Block_Sz), trunc_matrx)
        Block_Sz = trunc_matrx' * temp

        mul!(temp, kron(sparse(I, m, m), Block_Sp), trunc_matrx)
        Block_Sp = trunc_matrx' * temp

        mul!(temp, kron(sparse(I, m, m), Block_Sm), trunc_matrx)
        Block_Sm = trunc_matrx' * temp

        m = trunc_index # for next iterations
    end        
end

#Definition of starting parameters
Delta::Float64 = 0.5
max_states::Int16 = 5
NIter::Int32 = 5 #Final lattice size: 2^Iter

#output vectors
Energy = Vector{Float64}(undef, NIter)
Energy_2 = Vector{Float64}(undef, NIter)
Energy_bond = Vector{Float64}(undef, NIter)
Energy_bond[1] = 0 #no bond

NRG(Delta, max_states, NIter, Energy, Energy_2, Energy_bond)
for i = 1:NIter
    println(lpad(2^i, 2, "0"), " | ", Energy[i], " | ", Energy_bond[i], " | ", Energy_2[i])
end
1 Like
2

Try again with

evals, evecs = eigen!(Hermitian(H_dense_super))

It is likely that H_dense_super is not exactly hermitian and then the diagonalisation yields you complex numbers with very small imaginary parts (as the error suggests). If you wrap the matrix in Hermitian then hermiticity is ensured. (Effectively the lower triangular part (IIRC) is mirrored across the diagonal, so the matrix is exactly hermitian).

3

Well, you solved an issue, but not the one in this question, ahahhahahaha.
To wrap it up:
There were some little (non)zeros in the matrix and writing Symmetric() solved the issue.
Now numeric results go as expected.

However it solved it without using MKL.
Without MKL:

02 | -0.37499999999999944 | 0.0 | -0.18749999999999972
04 | -0.9098524898764083 | -0.26742624493820444 | -0.22746312246910208
08 | -1.9485400435360116 | -0.25967188841490085 | -0.24356750544200145
16 | -4.017274038426099 | -0.25859174936126095 | -0.2510796274016312
32 | -8.154573454879735 | -0.2585812135283522 | -0.2548304204649917

With MKL and MKLSparse:

02 | -0.37499999999999944 | 0.0 | -0.18749999999999972
04 | -0.8321067811865448 | -0.22855339059327268 | -0.2080266952966362
08 | -0.7678727796362532 | 0.016058500387572894 | -0.09598409745453165
16 | -0.09950923658373387 | 0.08354544288156492 | -0.006219327286483367
32 | -7.511857594535956e222 | -4.694910996584972e221 | -2.347455498292486e221

and while I know the first results are now right, the last ones are wrong

4

I think there are 2 errors in your question :wink:

I don’t know why MKL gives you wrong results though. Can you try just adding one of MKL and MKLSparse to better isolate the issue?

Also your error presents different now I think. The last line here is very wrong when before it was still of the same order of magnitude as the correct result.

5

Bug when using MKLSparse following use of MKL · Issue #36 · JuliaSparse/MKLSparse.jl · GitHub might be the reason.

See also Conflict between MKL_jll and libgomp in CompilerSupportLibraries_jll · Issue #700 · JuliaPackaging/BinaryBuilder.jl · GitHub.

1 Like
6

Do you still need an example from my program to pin-point the issue or should I wait this to be solved? I don’t know if I could make it simpler than the one linked and I surely don’t think I can make a general example, but just a specific one

7

I don’t think so.
I will write a difficult to read program because I know its analytical solution. It’s a little variation from the simpler one I wrote above. It’s a variational algorithm and it converges much faster.
The analitical result is -0.274519

Read the last result of the 4th column:

Without MKL and MKLSparse:
-2.699462703469799e-01
With MKL:
-2.699500319370028e-01 (a better one, but still: allarming number of changed digits)
With MKLSparse:
-2.614410653180197e-01 (what happened here?)
With both:
-1.534232626884432e-01 (?)

using MKL
using MKLSparse
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)

    Random.seed!(1)

    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

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

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("| %.15e | %.15e | %.15e | %.15e | \n", Energy[i], Energy_bond[i], Energy_2[i], Truncation_Error[i])
end