performance on parallels

[cfengno1]

Dear MPSKit Developer,

Thank you for this wonderful package!

I am trying to use multithreading in IDMR
Uploading Archive.zip…
G calculations, but it seems that the usage of cores is pretty low: only 3 out of 3 cores are used.

I have attached my scripts. Can you help to improve the parallelism?

Thank you!

Best,
Feng

[lkdvos]

Hi Feng,

It seems like you didn't upload the files here, and something went wrong.

Considering the parallelization in IDMRG, this is indeed not straightforward, and can depend on a variety of settings. Some of this information I can deduce once you upload the script, but mostly the questions are about what kind of model you are considering, and what kind of symmetry you are using.
For dense (non-symmetric) tensors, presumably the best thing to do is to hand most of the threads you have available to BLAS, which is just a setting of LinearAlgebra.BLAS itself, and possibly you could experiment if using threads for Strided.jl is beneficial as well.
For symmetric tensors the story changes slightly, and there are some hidden tensorkit features that could be worth playing around with, but it all depends quite a bit on the specific case.

Best,
Lukas

[cfengno1]

Dear Lukas,

Here is the code attached. Thank you so much!
Archive.zip

[cfengno1]

Dear Luka,

Could you help?

Best,
Feng

[lkdvos]

Sorry for being a little slow on this, I will definitely try and help but got sidetracked with some other projects, feel free to remind me again if I don't get to this by Friday, I'll do my best!

[cfengno1]

Thank you :)

[cfengno1]

Dear Lukas,

Sorry to bother you again. I hope you will have time to look into this. Thank you!

Best,
Feng

[lkdvos]

Hi Feng,

I finally got to looking at your code in a bit more detail, thanks for your patience!

Let me start with the long story short: the low core usage is not necessarily something you configured wrong, it is mostly a property of the IDMRG algorithm itself.
However, before getting into that, on a separate note there is some code in model_setup.jl that probably is not what you meant.

Hopping operator

In hopping_SU2 you build the creation operator as

hop_l = ones(elt, pspace ← pspace ⊗ aspace)    # c†
hop_r = ones(elt, aspace ⊗ pspace ← pspace)    # c

ones sets every reduced matrix element to 1.
With SU(2) spin symmetry those reduced matrix elements are not all 1: the Wigner–Eckart theorem fixes their relative values, and they come out as 1, -1, √2, √2.
Comparing your operator block-by-block against the correct one gives:

channel (codomain ← domain)	correct	yours	ratio
(1,½),(0,0) ← (0,0),(1,½)	1	1	1
(0,2),(0,0) ← (1,½),(1,½)	√2	1	0.707
(1,½),(1,½) ← (0,0),(0,2)	√2	1	0.707
(0,2),(1,½) ← (1,½),(0,2)	−1	−½	0.5

Since the ratios differ, this is not an overall rescaling of t — the doublon channels are weighted differently from the empty-site channel, so it is a genuinely different Hamiltonian.

The cleanest way to see the consequence is the free-electron limit. For a half-filled open chain of 20 sites at U = 0 the exact answer is E = 2 Σ_{k=1}^{10} (-2t cos(kπ/21)) = -24.76298. Running DMRG2 with both operators:

U = 0 (free electrons), L = 20, half filling
  exact (free fermions)        E = -24.76298000
  TensorKitTensors e_hopping   E = -24.76297821   (error 1.8e-06)
  your `ones` hopping          E = -17.86281469   (error 6.9e+00)

and at your actual U = 8, E/site = -0.3173 vs -0.1636.

I would strongly recommend not hand-writing these operators. TensorKitTensors.jl has them all, tested, for every symmetry combination:

using TensorKitTensors.HubbardOperators
P    = hubbard_space(U1Irrep, SU2Irrep)
hop  = e_hopping(ComplexF64, U1Irrep, SU2Irrep)   # already c†c + h.c.
d_op = ud_num(ComplexF64, U1Irrep, SU2Irrep)      # n↑ n↓
n_op = e_num(ComplexF64, U1Irrep, SU2Irrep)       # n↑ + n↓

Note e_hopping is already Hermitian, so it goes on each bond once — -t * hop{i,j}, not hop{i,j} + hop{j,i}.

One thing you did that is definitely worth keeping: shifting the U(1) charges to -1, 0, +1. With the natural labels N = 0, 1, 2 a half-filled infinite MPS cannot close its unit cell, since every cell would add +L charge. The supported way to express that is

I = sectortype(H)
H = MPSKit.add_physical_charge(H, fill(I(0, 1, 0), length(H)))

which is exactly the mechanism MPSKitModels' own n keyword uses.
That way you keep the standard operators and only relabel at the end.

2. Core usage

IDMRG2 sweeps sites strictly sequentially: the local eigensolve at site i+1 needs the environment produced at site i. There is no site-level parallelism to exploit, and correspondingly MPSKit's unit-cell multithreading (MPSKit.Defaults.set_scheduler!) is not used by IDMRG/IDMRG2 at all — it is only wired into VUMPS, VOMPS, TDVP, the Grassmann optimizer and the quasiparticle solver.

So julia -t 8 on its own buys you very little, and what you are seeing is probably just MKL_NUM_THREADS=4 — 4 cores inside one gemm at a time, which matches your "~3 out of 32" exactly.

For symmetric tensors the work is spread over many small blocks, one per symmetry sector, and in the current stack those loops are serial:

• mul! on a TensorMap loops over coupled sectors serially, and so do the block factorizations (foreachblock still has a # TODO: implement scheduler). Only BLAS runs inside each block.
• BlockTensorKit (which carries the MPO) has no threading.
• The permutations/braidings can be threaded over fusion blocks, but that is off by default.

That last one is the knob worth knowing about, which is somewhat hidden inside of TensorKit:

TensorKit.set_num_transformer_threads(n)     # threads permute/braid over fusion blocks
TensorKit.set_num_manipulation_threads(n)    # threads construction of recoupling data

Both default to 1, and both are capped by Threads.nthreads() — so these are the settings that actually make julia -t n do something in an IDMRG2 run. With SU(2) ⊠ U(1) ⊠ fℤ₂ and a 6-leg cylinder, permutes and recoupling are a substantial share of the time, so this might be where I would look first.

Concretely, I would scan the split of your 32 cores between the two pools rather than guess:

Nt=1  BLAS=32     # all cores to BLAS -- reasonable starting point
Nt=4  BLAS=8
Nt=8  BLAS=4      # what you have now

and set them explicitly in the script rather than through MKL_NUM_THREADS, because under Slurm Sys.CPU_THREADS reports the whole node and the automatic defaults will oversubscribe your allocation. Note also that MKL spawns its threads per Julia thread, unlike OpenBLAS — see the parallelism page — which is another reason not to crank both numbers at once.

Honestly though: at small bond dimension your blocks are simply too small for 32 cores, and no setting we currently have will fix that.

3. Other improvements

Your loop currently does

psi, _, _ = find_groundstate(psi, H, IDMRG2(...))
E = sum(real, expectation_value(psi, H)) / L

Both lines throw away work. find_groundstate returns the environments as its second output, and expectation_value(psi, H) without them solves the environment fixed-point equations from scratch — at χ = 8000 that is not cheap. Over 10 outer iterations you are paying for roughly 20 cold environment solves that you already had. Just thread them through:

envs = environments(psi, H, psi)
psi, envs, delta = find_groundstate(psi, H, IDMRG2(...), envs)
E = real(expectation_value(psi, H, envs)) / length(psi)

Passing envs back in also warm-starts each solve from the previous solution.

Ramp the bond dimension. Going straight to truncdim(8000) from a bond-dimension-1 start means the first ~10 sweeps are spent just growing χ geometrically, and then the sweeps that actually cost something get almost no iterations.
Something like 500 → 1000 → 2000 → 4000 → 8000 with more sweeps at the cheap stages converges much faster per unit wall time.

Checkpoint less often: Each find_groundstate call ends with a full re-gauge plus an environment recalculation, so maxiter = 2 pays that overhead every two sweeps. Chunks of 5–10 sweeps might be better trade. (IDMRG/IDMRG2 do not have the finalize callback that DMRG/VUMPS have, which is why you have to chunk at all — I will look into adding it.)
save(path * name, ...) in your script references path and name, which are never defined — the run would throw at the first checkpoint.

4. Change algorithm types

VUMPS, VOMPS and TDVP do parallelize over the unit cell through MPSKit.Defaults.set_scheduler!, because their local updates at different sites are independent. With a 6-site unit cell that is a real factor. A common pattern is to use IDMRG2 to grow the bond dimension and then switch to VUMPS at fixed χ to converge — worth trying if the IDMRG2 sweeps stay the bottleneck, or even the alternate between the two, using IDMRG2 as dynamical bond dimension but with VUMPS as the main convergence driver.

Best,
Lukas

[cfengno1]

Thank you Lukas! This is really helpful!