[BlockSparseArrays] BlockSparseArray functionality

This issue lists functionalities and feature requests for BlockSparseArray.

using LinearAlgebra
using NDTensors.BlockSparseArrays: BlockSparseArray

a = BlockSparseArray{Float64}([2, 3], [2, 3]);

Bugs

• LinearAlgebra.Adjoint{T,<:BlockSparseArray{T}} returns a BlockedArray in certain slicing operations (https://github.com/ITensor/BlockSparseArrays.jl/issues/2).
• Fix cat(::BlockSparseArray, ::BlockSparseArray) for dual axes (followup to ITensor/ITensors.jl#1579).

Feature requests

• Allow the syntax BlockSparseArray{Float64}([U1(0) => 2, U1(1) => 3], [U1(0) => 2, U1(1) => 3]) which could implicitly create axes with GradedUnitRange internally.
• Use NestedPermutedDimsArray instead of SparsePermutedDimsArrayBlocks (similar to how we are removing SparseAdjointBlocks/SparseTransposeBlocks in ITensor/ITensors.jl#1580). Started in ITensor/ITensors.jl#1589, ITensor/ITensors.jl#1590.
• An alternative design to ITensor/ITensors.jl#1589 would be to redefine NestedPermutedDimsArray as a PermutedDimsArray wrapping a MappedArray where the map and inverse map convert to PermutedDimsArray, that would be good to explore so we don't have to support all of the NestedPermutedDimsArrays code, which is mostly just a copy of Base.PermutedDimsArrays anyway.
• If slices are just blockwise, like b = @view a[Block.(1:2), [Block(2), Block(1)], define blocks(b) as @view blocks(a)[1:2, [2, 1]], as opposed to using the more general SparseSubArrayBlocks in those cases. Like the new NestedPermutedDimsArray, in principle SparseSubArrayBlocks could be replaced by a NestedSubArray type that defines the slicing behavior of the array storing the blocks and also the slicing of the blocks themselves, but that might be overkill and the concept is very particular to block arrays. But maybe SubArray of the blocks could still be used to simplify the code logic in SparseSubArrayBlocks.
• Constructor from Dictionary should check block sizes (https://github.com/ITensor/BlockSparseArrays.jl/issues/2).
• Blockwise linear algebra operations like svd, qr, etc. See ITensor/BlockSparseArrays.jl#3. These are well defined if the block sparse matrix has a block structure (i.e. the sparsity pattern of the sparse array of arrays blocks(a)) corresponding to a generalized permutation matrix. Probably they should be called something like block_svd, block_eigen, block_qr, etc. to distinguish that they are meant to be used on block sparse matrices with those structures (and error if they don't have that structure). See ^linalg for a prototype of a blockwise QR. See also BlockDiagonals.jl for an example in Julia of blockwise factorizations, they use a naming scheme svd_blockwise, eigen_blockwise, etc. The slicing operation introduced in ITensor/ITensors.jl#1489 will be useful for performing block-wise truncated factorizations.
• Define storedblockview(a, ::Block) to get the view of a stored block, which would allow avoiding wrapping the block in a view. Also define getstoredindex(a, ::Block).
• Change the behavior of slicing with non-blocked ranges (such as a[1:2, 1:2]) to output non-blocked arrays, and define @blocked a[1:2, 1:2] to explicitly preserve blocking. See the discussion in https://github.com/JuliaArrays/BlockArrays.jl/issues/347.
• Reconsider the design of how duality is stored in graded unit ranges (graded axes), for example storing it at the level of the sector labels with a new SectorDual type and/or as a boolean flag.
• Display non-initialized blocks differently from zero-initialized blocks, currently they both print as zeros.

Fixed

• Rename SparseArrayInterface to SparseArraysBase and move SparseArrayDOKs into the new SparseArraysBase. (Fixed in ITensor/ITensors.jl#1591, ITensor/ITensors.jl#1592.)
• Constructors like BlockSparseArray{Float64,2,Matrix{Float64}}([2, 3], [2, 3]), BlockSparseArray{Float64,2}([2, 3], [2, 3]), BlockSparseMatrix{Float64}([2, 3], [2, 3]) are not defined. (Fixed by ITensor/ITensors.jl#1586.)
• Rename block_nstored to block_stored_length and nstored to stored_length. (Fixed by ITensor/ITensors.jl#1585.)
• Rename BlockSparseArrayLike to AnyAbstractBlockSparseArray, which is the naming convention used in other Julia packages for a similar concept[^anygpuarray].
• Implement direct sums/concatenations of block sparse arrays that preserve and make use of block sparsity. This could be done by overloading Base.cat and related functions. (Implemented in ITensor/ITensors.jl#1579.)
• Fix issues with nested slices like @views a[[Block(2), Block(1)], [Block(2), Block(1)]][2:4, 2:4] in Julia 1.11 (see tests marked as broken in ITensor/ITensors.jl#1539). (Fixed in ITensor/ITensors.jl#1575.)
• Support for blocks that are on GPU. (Partially addressed in ITensor/ITensors.jl#1560 but more work is needed, we can track issues individually from now on.)
• Better support for blocks that are not Array, for example DiagonalArrays.DiagonalArray, SparseArrayDOKs.SparseArrayDOK, LinearAlgebra.Diagonal, etc. BlockSparseArray can have blocks that are AbstractArray subtypes, however some operations don't preserve those types properly (i.e. implicitly convert to Array blocks) or don't work. (Partially addressed in ITensor/ITensors.jl#1560 but more work is needed, we can track issues individually from now on.)
• Cannot create zero dimension array (BlockSparseArray{Float64}() fails). (Fixed in ITensor/ITensors.jl#1574.)
• permutedims crashes for some block sparse arrays (https://github.com/ITensor/BlockSparseArrays.jl/issues/2). (Fixed in ITensor/ITensors.jl#1574.)
• copy(adjoint) does not preserve dual axes. (Fixed in ITensor/ITensors.jl#1574.)
• block_stored_indices(::LinearAlgebra.Adjoint{T, BlockSparseArray}) does not transpose its indices. (Matt: I don't see an issue here, the values of block_stored_indices(a') are the nonzero/stored block locations, the keys of block_stored_indices(a') are an implementation detail and should not be used.)
• LinearAlgebra.norm(a) crashes when a contains NaN.
• Slicing a BlockSparseArray that has GradedUnitRange axes (as opposed to GradedOneTo) fails.
• Sub-slices of multiple blocks sometimes fails with a dual axis..
• a[:, :] creates an array with ill behaved axes (it should just be equivalent to copy(a)). Also triggers display error.
• Printing a BlockSparseArray triggers BoundsError in some rare contexts (https://github.com/ITensor/BlockSparseArrays.jl/issues/2).
• Automatically initializing a block when only a slice is used (this should work if you make sure to use a view, see https://github.com/ITensor/BlockSparseArrays.jl/issues/2).
• Implement slicing syntax for merging blocks, i.e. https://github.com/JuliaArrays/BlockArrays.jl/issues/359.
• a = BlockSparseArray{Float64}([2, 3], [2, 3]); @view a[Block(1, 1)] returns a SubArray where the last type parameter which marks whether or not the slice supports faster linear indexing is false, while it should be true if that is the case for that block of a (this is addressed by ITensor/ITensors.jl#1513, @view a[Block(1, 1)] no longer outputs a SubArray, but rather either the block data directly or a BlockView object if the block doesn't exist yet).
• TensorAlgebra.contract fails when called with view(::BlockSparseArray, ::Block) or reshape(view(::BlockSparseArray, ::Block), ...). As a workaround it will work if you use view!/@view! instroduced in ITensor/ITensors.jl#1498.
• Fix tests that are marked as broken in https://github.com/ITensor/ITensors.jl/blob/main/NDTensors/src/lib/BlockSparseArrays/test/test_basics.jl.
• a = BlockSparseArray{Float64}([2, 3], [2, 3]); b = @view a[Block.(1:2), Block.(1:2)]; b[Block(1, 1)] = randn(2, 2) doesn't set the block Block(1, 1) (it remains uninitialized, i.e. structurally zero). I think the issue is that @view b[Block(1, 1)] makes two layers of SubArray wrappers instead of flattening down to a single layer, and those two layers are not being dispatched on properly (in general we only catch if something is a BlockSparseArray or a BlockSparseArray wrapped in a single wrapper layer).
• Compatibility issues with BlockArrays v1.1, see CI for ITensor/ITensors.jl#1363. Fixed by ITensor/ITensors.jl#1503.
• Rewrite GradedAxes using BlockArrays v1[^blockarrays_v1].
• r = gradedrange([U1(0) => 1]); a = BlockSparseArray{Float64}(r, r); size(view(a, Block(1,1))[1:1,1:1]) returns a tuple of LabelledInteger instead of Int (see discussion, keep it that way at least for now).
• r = gradedrange([U1(0) => 1]); a = BlockSparseArray{Float64}(dual(r), r); @view(a[Block(1, 1)])[1:1, 1:1] and other combinations of dual lead to method ambiguity errors.
• Implement slicing syntax for slicing multiple subsets of blocks, i.e. https://github.com/JuliaArrays/BlockArrays.jl/issues/358.
• dual is not preserved when adding/subtracting BlockSparseArrays, i.e. g = gradedrange([U1(0) => 1]); m = BlockSparseArray{Float64}(dual(g), g); isdual(axes(m + m, 1)) should be true but is false.
• Error printing views of blocks of BlockSparseArray with GradedUnitRange axes, i.e.r = gradedrange([U1(0) => 1]); a = BlockSparseArray{Float64}(r, r); @view a[Block(1, 1)].
• Change the design of slicing with unit ranges, i.e. a[2:4, 2:4], by using BlockArrays.BlockSlice.
• Fix a[Block(2), Block(2)] = randn(3, 3).
• Fix a[Block(2, 2)] .= 1.
• Fix @view(a[Block(1, 1)])[1:1, 1:1] = 1.
• Throw error if a[Block(1, 1)] = b if size(a[Block(1, 1)]) != size(b).
• Fix matrix multiplication of BlockSparseMatrix involving dual axes.
• Fix matrix multiplication involving adjointed BlockSparseMatrix, i.e. a' * a and a * a', with and without dual axes.
• Take the dual of the axes in adjoint(::BlockSparseMatrix). Can be implemented by overloading axes(::Adjoint{<:Any,<:AbstractBlockSparseMatrix}).
• show(::Adjoint{<:Any,<:BlockSparseMatrix}) and show(::Transpose{<:Any,<:BlockSparseMatrix)) are broken.
• fix eachindex(::BlockSparseArray) involving dual axes.
• Fix [Hermitian] transposed BlockSparseMatrix, i.e. a' (in progress in^adjoint).
• Base.similar(a::BlockSparseArray, eltype::type) and Base.similar(a::BlockSparseArray, eltype::type, size::NTuple{N,AbstractUnitRange}) do not set eltype
• Make copy(::BlockSparseArray) copy the blocks.
• Slicing with a[1:2, 1:2] is not implemented yet and needs to be implemented (in progress in^slicing).
• Function for accessing a list of initialized blocks. Currently you can use stored_indices(blocks(a))) to get a list of Block corresponding to initialized/stored blocks. Ideally there would be shorthands for this like block_stored_indices(a) (in progress in^slicing).
• Function for accessing the number of initialized blocks. Currently you can use nstored(blocks(a)) to get the number of initialized/stored blocks. Ideally there would be shorthands for this like block_nstored(a) (in progress in^slicing).
• In place operations .*= and ./=, such asa .*= 2, are broken (in progress in[^1]).
• Base.:*(::BlockSparseArray, x::Number) and Base.:/(::BlockSparseArray, x::Number) are not defined
• Base.:*(::ComplexF64, ::BlockSparseArray{Float64}) does not change data type for empty array and crashes if a contains data.

[^blockarrays_v1]: https://github.com/ITensor/ITensors.jl/pull/1452, https://github.com/JuliaArrays/BlockArrays.jl/pull/255
[^anygpuarray]: https://github.com/JuliaGPU/GPUArrays.jl/blob/v11.1.0/lib/GPUArraysCore/src/GPUArraysCore.jl#L27, https://github.com/JuliaGPU/CUDA.jl/blob/v5.4.2/src/array.jl#L396

This is a broad tracker for BlockSparseArray, covering multiple unchecked bugs and feature requests rather than one scoped change. Start by selecting a single unchecked item, then read the linked issue or reference, such as NDTensors/src/lib/BlockSparseArrays/test/test_basics.jl where mentioned. Done means the selected behavior is implemented and its relevant tests pass.