Hello, thanks for replying

since the example uses nearest neighbor gates, I stripped away all the extra logic from the apply function and wrote this

```
function applyNNgates(gates::Vector{ITensor},psi::ITensors.AbstractMPS; kwargs...)
orthogonalize!(psi, 1)
for gate in gates
# apply gate
ns=findsites(psi,gate)
ns = sort(ns)
wf = (psi[ns[1]] * psi[ns[2]]) * gate
noprime!(wf)
# restore MPS form
linds = uniqueinds(psi[ns[1]],psi[ns[2]])
# L,R = factorize(wf,linds;ortho="left",kwargs...) <- uses eigen(...), bare implementation below
simlinds = sim(linds)
wf2 = wf * replaceinds(dag(wf), linds, simlinds)
F = eigen(wf2, linds, simlinds; ishermitian=true, kwargs...)
L = F.Vt
R = dag(L) * wf
psi[ns[1]] = L
psi[ns[2]] = R
# move ortho center
orthogonalize!(psi, ns[1])
end
return psi
end
```

I still get the same slow execution time. The extra logic from the apply() and factorize() functions doesn’t seem to have much effect. It is mostly if else statements and no swapping is happening in the example.

For the system (Heisenberg, N=100, tau=0.001 and ttotal=1) I timed the calls on both languages and found this discrepancy

Julia: orthogonalize() ~ 0.2 ms

c++ : psi.position() ~ 0.04 ms

orthogonalize() is a major source of slowdown. I checked with a profiler if the boilerplate around svd() is causing the slowdown, the answer is no. I also changed the algorithm to “recursive” which gave a slightly slower time.

# About using larger MPS bond dimensions

I ran into a different problem while attempting this. I’m again using the original time evolution code from Julia and C++ documentation.

let’s take for example N=32, cutoff=1E-32, tau = 0.001, ttotal = 1

I monitor the MPS bond dimension on both languages and something strange happens; the bond dimensions is exactly equal up until t=0.625

Julia:

…

(maxlinkdim(psi), t) = (40, 0.622)

(maxlinkdim(psi), t) = (40, 0.623)

(maxlinkdim(psi), t) = (53, 0.624)

(maxlinkdim(psi), t) = (74, 0.625) <<<<

(maxlinkdim(psi), t) = (92, 0.626)

(maxlinkdim(psi), t) = (106, 0.627)

(maxlinkdim(psi), t) = (119, 0.628)

(maxlinkdim(psi), t) = (124, 0.629)

(maxlinkdim(psi), t) = (129, 0.63)

…

// finishes time evolution ~ 3 mins

C++:

…

maxbond= 40 t: 0.622

maxbond= 40 t: 0.623

maxbond= 53 t: 0.624

maxbond= 130 t: 0.625 <<<<

maxbond= 288 t: 0.626

maxbond= 666 t: 0.627

// freezes

Also, I’m not sure if a bond dimension of ~40 is enough to drown out any overhead, but it’s worth mentioning that c++ was much faster than julia up until the freeze.