LinearMap Error

LinearMapErr.jl

using LinearMaps, LinearAlgebra
"""
Ising model
"""
function buildIsing(h=1.0) # Ising model with transverse magnetic field h (critical h=1 by default)
    X = [0 1; 1 0]
    Z = [1 0; 0 -1]
    Iden = Matrix(1.0I, (2, 2))
    XX = kron(X, X)
    ZI = kron(Z, Iden)
    IZ = kron(Iden, Z)
    H2 = -(XX + h/2 * (ZI + IZ))
    return H2
end

"""
Multiply local hamiltonian h with the total wavefunction.
See Vidal's explanation in PSI.
"""
function multiplyHPsi!(nPsi::AbstractVector, Psi::AbstractVector)  # requires H2 to have been defined
    h = 1
    H2 = buildIsing(h) # Hamiltonian
    D,U = eigen(H2)
    shiftE = D[end]
    # remember to re-shift the energy later on!
    H2 = H2 - shiftE * Matrix(1.0I, (4, 4))
    # input and output wavefunction should have the same dimensions
    length(nPsi)==length(Psi) || throw(DimensionMismatch())
    dimN = length(Psi)  # Dimension of the vector space
    N = convert(Int64, log2(dimN))  # Number of spins
    dimPsi = tuple(fill(2::Int,(1, N))...) # dimsPsi =(2,2,2, ..., 2)
    # for shifting the spin position by 1
    a = collect(2:N)
    a = [a;1]
    # tuple for cyclic permutation cyclperm = (2, 3, ... N, 1)
    cyclperm = tuple(a...) 
    Psi = reshape(Psi, dimPsi)
    nPsi = zeros(dimPsi)
    for n in 1:N
        # multplication of H2 using matrix shape
	Psi = reshape(Psi, (4, 2^(N-2)))
	nPsi = reshape(nPsi, (4, 2^(N-2)))
	# Act on H2 to the input wavefunction and save
	# it to the output wavefunction
	nPsi += H2 * Psi
	# shift the leg using the tensor shape
	Psi = reshape(Psi, dimPsi)
        nPsi = reshape(nPsi, dimPsi)
        Psi = permutedims(Psi, cyclperm)
        nPsi = permutedims(nPsi, cyclperm)
    end
    nPsi = reshape(nPsi, dimN)
    return nPsi
end

# Define the LinearMap from `multiplyHPsi!`
N = 2  # number of spins
H = LinearMap{Float64}(multiplyHPsi!, 2^N; ismutating=true, issymmetric=true)
# test the actions of `H` and `multiplyHPsi!` on the same random vector, they should give the same result
# initialize the input and output vector `Psi` and `nPsi`
Psi = randn(2^N)
nPsi = randn(2^N)
# normalize the input vector
Psi = Psi / norm(Psi)
# action of `multiplyHPsi!`. It looks correct
nPsi1 = multiplyHPsi!(nPsi, Psi)
# action of `H`. It doesn't look correct
nPsi2 = H(Psi)
println("Output vector after the action of multiplyHPsi! is:")
println(nPsi1)
println("Output vector after the action of H is:")
println(nPsi2)