Comparison of eig solution between MATLAB and Julia for the TFIM

eigbench.jl

using ProfileView, Gadfly

function TFIM_1D_PBC(N,J,H)
# Hilbert Space Dimension
dim = 2^N
# Hamiltonian Matrix
ham = zeros(dim,dim)
# Build the diagonal part of the Hamiltonian
# First loop through the hilbert space
for bra in 0:dim-1
  # Next, loop through all the Hamiltonian elements
  diag_term = 0
  for s in 0:N-1 # Take site "s" (up to the last spin)
    S0 = 2*((bra>>s)&1) - 1 # Look at spin of site "s"
    S1 = 2*((bra>>((s+1)%N))&1) - 1 # Look at spin of site "s+1 mod N"
    diag_term = diag_term - J*S0*S1
  end
  ham[bra+1,bra+1] = diag_term
  # Now let's do the off-diagonal terms
  for s in 0:N-1
    y = 2^s # This labels the bit to be flipped
    ket = bra $ y # Flip that bit in the bra
    ham[bra+1,ket+1] = -H
  end
end
return ham
end

function diagonalize()
N = 8; J = 1
hrange = [0:0.001:0.1;0.2:0.1:2] # More points for lower h/J ratios
E = zeros(length(hrange),2^N)
V = zeros(length(hrange),2^N,2^N)
for (i,h) in enumerate(hrange)
    H = TFIM_1D_PBC(N,J,h)
    a,b = eig(H) # Diagonalize
    if (h<0.03)&&(abs(b[1,1])<abs(b[2^N,1]))
        # Track the correct energy level since numerical errors sometime
        # cause 1 & 2 to be flipped)
        temp = b[:,1]
        b[:,1] = b[:,2]
        b[:,2] = temp
    end
    E[i,:] = a/N # Store the energy per spin for every eigenstate
    V[i,:,:] = b # Store the eigenvectors
end
    return E,V
end

E,V = @profile diagonalize()
ProfileView.view()

p = plot(x=1:50,y=abs(V[1:50,1,1]),Geom.line)
draw(PNG("eigbench_jl.png", 12cm, 10cm), p)

eigbench.m

function eigbench( )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
tic
[hrange,E,V] = diagonalize();
toc

plot(abs(V(1:50,1,1)))
export_fig('eigbench_mat.png')
end

function [ hrange,E,V ]= diagonalize()
N = 8; J = 1;
hrange = [0:0.001:0.1,0.2:0.1:2]; % More points for lower h/J ratios
E = zeros(length(hrange),2^N);
V = zeros(length(hrange),2^N,2^N);
for i = 1:length(hrange)
    h = hrange(i);
    H = TFIM_1D_PBC(N,J,h);
    [b,a] = eig(H); % Diagonalize
    if (h<0.03)&&(abs(b(1,1))<abs(b(2^N,1)))
        % Track the correct energy level since numerical errors sometime
        % cause 1 & 2 to be flipped)
        temp = b(:,1);
        b(:,1) = b(:,2);
        b(:,2) = temp;
    end
    E(i,:) = diag(a)/N; % Store the energy per spin for every eigenstate
    V(i,:,:) = b; % Store the eigenvectors
end
end

function ham = TFIM_1D_PBC(N,J,H)
% Hilbert Space Dimension
dim = 2^N;
% Hamiltonian Matrix
ham = zeros(dim,dim);
%Build the diagonal part of the Hamiltonian
%First loop through the hilbert space
for bra = 0:dim-1
    % Next, loop through all the Hamiltonian elements
    diag_term = 0;
    for s = 0:N-1 % Take site "s" (up to the last spin)
        S0 = 2*(mod(bitshift(bra,-s),2)==1) - 1; % Look at spin of site "s"
        S1 = 2*(mod(bitshift(bra,-mod(s+1,N)),2)==1) - 1; % Look at spin of site "s+1 mod N"
        diag_term = diag_term - J*S0*S1;
    end
    ham(bra+1,bra+1) = diag_term;
    %Now let's do the off-diagonal terms
    for s = 0:N-1
        y = 2^s;  % This labels the bit to be flipped
        ket = bitxor(bra,y);     %Flip that bit in the bra
        ham(bra+1,ket+1) = -H;
    end
end
end