[BlockSparseArrays] Function for fusing dimensions of BlockSparseArray

[mtfishman]

This introduces fusedims, a function for fusing arbitrary groups of dimensions of a BlockSparseArray. This is a superset of combiner contraction, since it allows combining multiple groups of indices in one operation (for example, it can matricize a block sparse array in one operation, i.e. for use in ITensor SVD code).

This isn't being used inside combiner contraction yet but the goal will be to use this as a backend for that.

It additionally adds some other useful block sparse operations like slicing according to subblocks. Here is an example of both functionalities:

using NDTensors.BlockSparseArrays: BlockSparseArray, BlockRange, Block, gradedrange, fusedims
using ITensors: QN

d = [1, 1]
sectors = [QN(0, 2), QN(1, 2)]
i = gradedrange(d, sectors)

B = BlockSparseArray{Float64}(i, i, i, i)
for b in BlockRange(B)
  if iseven(sum(b.n))
    B[b] = randn(1, 1, 1, 1)
  end
end

# Block slicing
B_sub = B[[Block(2)], [Block(2)], [Block(2)], [Block(2)]]

# Fusion
B_fused = fusedims(B, (1, 2), (3, 4))

This outputs:

julia> B
2×2×2×2-blocked 2×2×2×2 BlockSparseArray{Float64, 4, Array{Float64, 4}, NDTensors.BlockSparseArrays.SparseArray{Array{Float64, 4}, 4, NDTensors.BlockSparseArrays.BlockZero{NTuple{4, NDTensors.BlockSparseArrays.GradedBlockedUnitRange{Vector{Int64}, QN}}}}, NTuple{4, NDTensors.BlockSparseArrays.GradedBlockedUnitRange{Vector{Int64}, QN}}}:
[:, :, 1, 1] =
 -0.359803   0.0
  0.0       -0.757009

[:, :, 2, 1] =
  0.0      -0.468575
 -0.87562   0.0

[:, :, 1, 2] =
 0.0     -0.592078
 1.8012   0.0

[:, :, 2, 2] =
 -1.80674   0.0
  0.0      -0.444399

julia> B_sub
1×1×1×1-blocked 1×1×1×1 BlockSparseArray{Float64, 4, Array{Float64, 4}, NDTensors.BlockSparseArrays.SparseArray{Array{Float64, 4}, 4, NDTensors.BlockSparseArrays.BlockZero{NTuple{4, BlockArrays.BlockedUnitRange{Vector{Int64}}}}}, NTuple{4, BlockArrays.BlockedUnitRange{Vector{Int64}}}}:
[:, :, 1, 1] =
 -0.4443991343705862

julia> B_fused
2×2-blocked 4×4 BlockSparseArray{Float64, 2, Matrix{Float64}, NDTensors.BlockSparseArrays.SparseArray{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}:
 -0.359803  -1.80674   │   0.0        0.0     
 -0.757009  -0.444399  │   0.0        0.0     
 ──────────────────────┼──────────────────────
  0.0        0.0       │  -0.87562    1.8012  
  0.0        0.0       │  -0.468575  -0.592078

@emstoudenmire

[mtfishman]

Note that the BlockSparseArray code doesn't know anything about the ITensors.QN type, it has a generic type for associating data with blocks on each dimension (called GradedBlockedUnitRange, with a constructor function gradedrange) and uses generic code to fuse QNs (assuming they can be added), sort them, and therefore determine the block permutations and mergers needed to perform abelian fusion.

[mtfishman]

There are a few limitations right now:

1. Currently the output of fusedims has all blocks allocated (instead of just the nonzero ones), that would have a simple fix but I just wanted to get this merged in a working state and make adjustments like that later.
2. fusedims drops sector information, again that is a simple fix which I will do later (it just requires making more convenient BlockSparseArray constructors to use in fusedims).
3. I haven't implemented uncombining but in principle that can use a lot of the same code logic.
4. There are some other optimizations we could do like doing the tensor permutation and slicing lazily, which is effectively what the current combiner contraction code in NDTensors does which is one reason why it is fairly complicated code. However, having the AbstractArray/BlockArray abstraction layer makes it much easier to handle the permutation lazily, so in principle it should be a small adjustment of the code introduced in this PR (i.e. replace permutedims with PermutedDimsArray and overload a few operations on PermutedDimsArray for operations like block slicing or copying as needed).