Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[BlockSparseArrays] More general broadcasting and slicing #1332

Merged
merged 19 commits into from
Mar 22, 2024

Conversation

mtfishman
Copy link
Member

@mtfishman mtfishman commented Feb 9, 2024

This is work in progress for defining map! for BlockSparseArrays with mismatched blocking.

This is functionality that is required to perform block fusion, as a follow-up to #1326.

So far, it performs some unnecessary copying, and also only works with one input array, which both should be relatively easy fixes.

  • Add tests.
  • For out-of-place expressions like a + b, currently it outputs a block sparse array with axes(a). If the axes/blocking of a and b are different, it should combine the axes (i.e. output a block sparse array with a finer blocking structure compatible with the blocking of a and b).
  • map(+, a, b) outputs a dense block array, not a block sparse array.
  • Fix issues with using GradedAxes as axes of block sparse arrays and performing mapping/broadcasting operations. There are some lingering bugs with that, and we need to make sure we preserve the sectors of the graded axes properly.
  • Generalize new map! functionality to arbitrary number of input arrays.
  • Avoid extraneous copying of source blocks (requires fixing a bug when using views of a BlockSparseArray).
  • Investigate some other slicing errors I came across during this PR, either fix or add broken tests.

A demonstration of new functionality:

using BlockArrays: Block, blockedrange, blocksize
using NDTensors.BlockSparseArrays: BlockSparseArray
using Random: randn!

function blockdiagonal(f!, elt::Type, axes::Tuple)
  a = BlockSparseArray{elt}(axes)
  for i in 1:minimum(blocksize(a))
    b = Block(ntuple(Returns(i), ndims(a)))
    a[b] = f!(a[b])
  end
  return a
end

elt = Float64
d_src = blockedrange([4, 4])
d_dest = blockedrange([2, 4, 2])
a_src = blockdiagonal(randn!, elt, (d_src, d_src))
a_dest = blockdiagonal(randn!, elt, (d_dest, d_dest))

a_dest .= 10 .* a_src

which outputs:

julia> a_src
2×2-blocked 8×8 BlockSparseArray{Float64, 2, Matrix{Float64}, NDTensors.SparseArrayDOKs.SparseArrayDOK{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.405138   0.627248   1.26039     0.1523310.0        0.0        0.0        0.0     
  1.26651   -1.318      0.90538    -0.6692140.0        0.0        0.0        0.0     
 -0.231649  -1.14072   -0.0183986   0.6414010.0        0.0        0.0        0.0     
  0.398392   1.15866   -0.840044   -0.7749310.0        0.0        0.0        0.0     
 ─────────────────────────────────────────────┼────────────────────────────────────────────
  0.0        0.0        0.0         0.0-0.149512  -0.415454  -0.337988  -0.214309
  0.0        0.0        0.0         0.0-1.08143   -0.88407    2.10251    0.180943
  0.0        0.0        0.0         0.0-0.265952  -0.441355   1.02287   -0.343253
  0.0        0.0        0.0         0.0-1.94984    0.913825  -0.61027   -0.232461

julia> a_dest
3×3-blocked 8×8 BlockSparseArray{Float64, 2, Matrix{Float64}, NDTensors.SparseArrayDOKs.SparseArrayDOK{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}:
  1.89965   -0.9570950.0        0.0        0.0        0.00.0       0.0      
 -0.366502   0.09682340.0        0.0        0.0        0.00.0       0.0      
 ───────────────────────┼───────────────────────────────────────────────┼──────────────────────
  0.0        0.00.366896   0.337547  -0.698596  -0.2788810.0       0.0      
  0.0        0.0-0.128684  -0.14473    0.751413  -0.5615580.0       0.0      
  0.0        0.0-0.369947   0.229541   0.949884  -0.9966510.0       0.0      
  0.0        0.0-0.422792  -0.339941   0.322883  -0.06763050.0       0.0      
 ───────────────────────┼───────────────────────────────────────────────┼──────────────────────
  0.0        0.00.0        0.0        0.0        0.00.546015  0.0357726
  0.0        0.00.0        0.0        0.0        0.0-0.763948  1.39853  

julia> a_dest .= 10 .* a_src
3×3-blocked 8×8 BlockSparseArray{Float64, 2, Matrix{Float64}, NDTensors.SparseArrayDOKs.SparseArrayDOK{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}:
 -4.05138    6.2724812.6039     1.52331     0.0        0.00.0       0.0    
 12.6651   -13.189.0538    -6.69214     0.0        0.00.0       0.0    
 ─────────────────────┼───────────────────────────────────────────────┼────────────────────
 -2.31649  -11.4072-0.183986   6.41401    -0.698596  -0.2788810.0       0.0    
  3.98392   11.5866-8.40044   -7.74931     0.751413  -0.5615580.0       0.0    
  0.0        0.0-0.369947   0.229541   -1.49512   -4.15454-3.37988  -2.14309
  0.0        0.0-0.422792  -0.339941  -10.8143    -8.840721.0251    1.80943
 ─────────────────────┼───────────────────────────────────────────────┼────────────────────
  0.0        0.00.0        0.0        -2.65952   -4.4135510.2287   -3.43253
  0.0        0.00.0        0.0       -19.4984     9.13825-6.1027   -2.32461

so it can handle cases where the blocking doesn't line up, while before this would have errored.

It works by first creating axes with a finer blocking structure that lines up with the blocking structure of the destination and source arrays, and then applies block-wise operations using that finer blocking structure.

@codecov-commenter
Copy link

codecov-commenter commented Feb 18, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 53.78%. Comparing base (f4ad958) to head (38c717f).
Report is 1 commits behind head on main.

❗ Current head 38c717f differs from pull request most recent head bf5e92c. Consider uploading reports for the commit bf5e92c to get more accurate results

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files
@@             Coverage Diff             @@
##             main    ITensor/ITensors.jl#1332       +/-   ##
===========================================
- Coverage   84.40%   53.78%   -30.62%     
===========================================
  Files         100       99        -1     
  Lines        8581     8528       -53     
===========================================
- Hits         7243     4587     -2656     
- Misses       1338     3941     +2603     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@mtfishman mtfishman changed the title [WIP][BlockSparseArrays] map! with mismatched blocking [BlockSparseArrays] map! with mismatched blocking Mar 22, 2024
@mtfishman mtfishman changed the title [BlockSparseArrays] map! with mismatched blocking [BlockSparseArrays] More general map! Mar 22, 2024
@mtfishman mtfishman marked this pull request as ready for review March 22, 2024 01:58
@emstoudenmire
Copy link
Collaborator

Very cool!

@mtfishman
Copy link
Member Author

@ogauthe @emstoudenmire I'm merging this.

The status is that you should be able to perform arbitrary broadcasting operations (adding, scalar multiplication, permutations) of block sparse arrays, as well as perform pretty general slicing operations, and it should preserve block structures and block labels stored on the axes (including symmetry labels). It also needs more tests, but I want to merge this now as a starting point for future work. I'm sure there will be corner cases and bugs to work out as it is put into practice but I think it is a good starting point.

@mtfishman mtfishman changed the title [BlockSparseArrays] More general map! [BlockSparseArrays] More general broadcasting and slicing Mar 22, 2024
@mtfishman mtfishman merged commit 957f2af into main Mar 22, 2024
9 checks passed
@mtfishman mtfishman deleted the BlockSparseArrays_map_mismatched_blocking branch March 22, 2024 18:55
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants