Optimizations for block sparse operations

[mtfishman]

Here is an issue listing potential improvements to block sparse tensor operations.

• Add multithreading to contraction and combining.
• Use a Dict instead of a sorted Vector for the block offsets list (like Store block-offset list as an std::unordered_map ITensor#336).
• In contraction, try combining indices first and then contracting with a matrix-matrix multiplication.
• Try Jutho's Strided.jl package for faster permutations.

[mtfishman]

I was curious about how the strategy of combining and then contracting would work right now. Here is the function for testing that:

function contract_combine(A::ITensor, B::ITensor)
  A_uniqueinds = uniqueinds(A, B)
  B_uniqueinds = uniqueinds(B, A)
  AB_commoninds = commoninds(A, B)
  CA = combiner(A_uniqueinds)
  CB = combiner(B_uniqueinds)
  CAB = combiner(AB_commoninds)
  AC = A * CA * CAB
  BC = B * CB * dag(CAB)
  CC = AC * BC
  C = CC * dag(CA) * dag(CB)
  return C
end

Running the following code:

function main()
  Nmax = 8
  for N in 1:Nmax
    println("order = ", 2 * N)
    d = 1
    i = Index(QN(0,2) => d, QN(1,2) => d)
    is = IndexSet(ntuple(n -> settags(i, "i$n"), Val(N)))
    A = randomITensor(is'..., dag(is)...)
    B = randomITensor(is'..., dag(is)...)
    println("Contract:              ", @elapsed A' * B)
    println("Combine then contract: ", @elapsed contract_combine(A', B))
    println()
  end
end

I get:

order = 2
Contract:              7.1355e-5
Combine then contract: 0.004867469

order = 4
Contract:              0.000149889
Combine then contract: 0.006314275

order = 6
Contract:              0.000469303
Combine then contract: 0.006801207

order = 8
Contract:              0.00263369
Combine then contract: 0.007549091

order = 10
Contract:              0.027489216
Combine then contract: 0.01168987

order = 12
Contract:              0.525102576
Combine then contract: 0.034727689

order = 14
Contract:              4.987328203
Combine then contract: 0.146389961

order = 16
Contract:              49.175892434
Combine then contract: 2.307021707

where "order = N" means two order-N ITensors are contracted over N/2 shared indices. You can see there is some overhead to combining, and it currently only pays off when contracting order-10 ITensors with each other. The QN combiner code is not particularly optimized right now, so I believe that could be improved a good amount. Clearly the scaling of not combining gets quite bad as the order increases (since there is an explosion in the number of nonzero blocks, I believe 2^(N-1) nonzero blocks, while after combining there are only 2 blocks).

[emstoudenmire]

Very interesting! Would be neat to see if later on it could become competitive even at order 6 or 4