From d3731434e7ded6d7622f67510442fda1dd287653 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Mon, 2 Jun 2025 15:59:51 -0400 Subject: [PATCH 01/11] Fix accidental method overwrite of default_svd_algorithm --- Project.toml | 2 +- src/factorizations/svd.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index c15f9724..327cea0e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BlockSparseArrays" uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" authors = ["ITensor developers and contributors"] -version = "0.6.8" +version = "0.6.9" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl index 187795ea..c73b8e24 100644 --- a/src/factorizations/svd.jl +++ b/src/factorizations/svd.jl @@ -12,7 +12,7 @@ struct BlockPermutedDiagonalAlgorithm{A<:MatrixAlgebraKit.AbstractAlgorithm} <: alg::A end -function MatrixAlgebraKit.default_svd_algorithm(A; kwargs...) +function MatrixAlgebraKit.default_svd_algorithm(A::AbstractBlockSparseMatrix; kwargs...) blocktype(A) <: StridedMatrix{<:LinearAlgebra.BLAS.BlasFloat} || error("unsupported type: $(blocktype(A))") # TODO: this is a hardcoded for now to get around this function not being defined in the From 37c843f90c57b5c383953b00ac27475d6a10fef9 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 3 Jun 2025 12:24:08 -0400 Subject: [PATCH 02/11] [WIP] Logical indexing --- Project.toml | 2 +- .../BlockArraysExtensions.jl | 19 ++++++++-- src/BlockArraysExtensions/blockedunitrange.jl | 35 +++++++++++++++++++ .../unblockedsubarray.jl | 5 ++- .../wrappedabstractblocksparsearray.jl | 6 ++++ .../blocksparsearrayinterface.jl | 8 +++++ 6 files changed, 71 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 327cea0e..ef35add8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BlockSparseArrays" uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" authors = ["ITensor developers and contributors"] -version = "0.6.9" +version = "0.6.10" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/src/BlockArraysExtensions/BlockArraysExtensions.jl b/src/BlockArraysExtensions/BlockArraysExtensions.jl index b4c2d926..3f811347 100644 --- a/src/BlockArraysExtensions/BlockArraysExtensions.jl +++ b/src/BlockArraysExtensions/BlockArraysExtensions.jl @@ -67,12 +67,19 @@ for f in (:axes, :unsafe_indices, :axes1, :first, :last, :size, :length, :unsafe @eval Base.$f(S::BlockIndices) = Base.$f(S.indices) end Base.getindex(S::BlockIndices, i::Integer) = getindex(S.indices, i) + +function _blockslice(x, y::AbstractUnitRange{<:Integer}) + return BlockSlice(x, y) +end +function _blockslice(x, y::AbstractVector{<:Integer}) + return BlockIndices(x, y) +end function Base.getindex(S::BlockIndices, i::BlockSlice{<:Block{1}}) # TODO: Check that `i.indices` is consistent with `S.indices`. # It seems like this isn't handling the case where `i` is a # subslice of a block correctly (i.e. it ignores `i.indices`). @assert length(S.indices[Block(i)]) == length(i.indices) - return BlockSlice(S.blocks[Int(Block(i))], S.indices[Block(i)]) + return _blockslice(S.blocks[Int(Block(i))], S.indices[Block(i)]) end # This is used in slicing like: @@ -347,7 +354,7 @@ function blockrange( axis::AbstractUnitRange, r::BlockVector{<:BlockIndex{1},<:AbstractVector{<:BlockIndexRange{1}}}, ) - return map(b -> Block(b), blocks(r)) + return map(Block, blocks(r)) end # This handles slicing with `:`/`Colon()`. @@ -365,6 +372,14 @@ function blockrange(axis::AbstractUnitRange, r::AbstractVector{<:Integer}) return Block.(Base.OneTo(1)) end +function blockrange(axis::AbstractUnitRange, r::BlockIndexVector) + return Block(r):Block(r) +end + +function blockrange(axis::AbstractUnitRange, r::Vector{<:BlockIndexVector}) + return map(Block, r) +end + function blockrange(axis::AbstractUnitRange, r) return error("Slicing not implemented for range of type `$(typeof(r))`.") end diff --git a/src/BlockArraysExtensions/blockedunitrange.jl b/src/BlockArraysExtensions/blockedunitrange.jl index f685b761..eb05e4dc 100644 --- a/src/BlockArraysExtensions/blockedunitrange.jl +++ b/src/BlockArraysExtensions/blockedunitrange.jl @@ -10,6 +10,7 @@ using BlockArrays: BlockVector, block, blockedrange, + blockfirsts, blockindex, blocklengths, findblock, @@ -169,6 +170,19 @@ function blockedunitrange_getindices( return mortar(map(b -> a[b], blocks(indices))) end +function blockedunitrange_getindices( + a::AbstractBlockedUnitRange, indices::AbstractArray{Bool} +) + blocked_indices = BlockedVector(indices, axes(a)) + return mortar( + map(Base.OneTo(blocklength(blocked_indices))) do b + binds = blocked_indices[Block(b)] + bstart = blockfirsts(only(axes(blocked_indices)))[b] + return findall(binds) .+ (bstart - 1) + end, + ) +end + # TODO: Move this to a `BlockArraysExtensions` library. function blockedunitrange_getindices(a::AbstractBlockedUnitRange, indices) return error("Not implemented.") @@ -197,6 +211,27 @@ function to_blockindices(a::AbstractBlockedUnitRange{<:Integer}, I::UnitRange{<: ) end +struct BlockIndexVector{T<:Integer,I<:AbstractVector{T},TB<:Integer} <: + AbstractVector{BlockIndex{1,Tuple{TB},Tuple{T}}} + block::Block{1,TB} + indices::I +end +Base.length(a::BlockIndexVector) = length(a.indices) +Base.size(a::BlockIndexVector) = (length(a),) +BlockArrays.Block(a::BlockIndexVector) = a.block +Base.getindex(a::BlockIndexVector, I::Integer) = Block(a)[a.indices[I]] +Base.copy(a::BlockIndexVector) = BlockIndexVector(a.block, copy(a.indices)) + +function to_blockindices(a::AbstractBlockedUnitRange{<:Integer}, I::AbstractArray{Bool}) + I_blocks = blocks(BlockedVector(I, blocklengths(a))) + return mortar( + map(eachindex(I_blocks)) do b + I_b = findall(I_blocks[b]) + BlockIndexVector(Block(b), I_b) + end, + ) +end + # This handles non-blocked slices. # For example: # a = BlockSparseArray{Float64}([2, 2, 2, 2]) diff --git a/src/abstractblocksparsearray/unblockedsubarray.jl b/src/abstractblocksparsearray/unblockedsubarray.jl index 355ae824..fc80e92f 100644 --- a/src/abstractblocksparsearray/unblockedsubarray.jl +++ b/src/abstractblocksparsearray/unblockedsubarray.jl @@ -4,7 +4,10 @@ using BlockArrays: BlockArrays, Block, BlockIndexRange, BlockSlice using TypeParameterAccessors: TypeParameterAccessors, parenttype, similartype const UnblockedIndices = Union{ - Vector{<:Integer},BlockSlice{<:Block{1}},BlockSlice{<:BlockIndexRange{1}} + Vector{<:Integer}, + BlockSlice{<:Block{1}}, + BlockSlice{<:BlockIndexRange{1}}, + BlockSlice{<:BlockIndexVector}, } const UnblockedSubArray{T,N} = SubArray{ diff --git a/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl b/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl index 84c24ed8..91705772 100644 --- a/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl +++ b/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl @@ -40,6 +40,12 @@ function Base.to_indices( return @interface BlockSparseArrayInterface() to_indices(a, inds, I) end +function Base.to_indices( + a::AnyAbstractBlockSparseArray, inds, I::Tuple{AbstractArray{Bool},Vararg{Any}} +) + return @interface BlockSparseArrayInterface() to_indices(a, inds, I) +end + # a[[Block(2), Block(1)], [Block(2), Block(1)]] function Base.to_indices( a::AnyAbstractBlockSparseArray, inds, I::Tuple{Vector{<:Block{1}},Vararg{Any}} diff --git a/src/blocksparsearrayinterface/blocksparsearrayinterface.jl b/src/blocksparsearrayinterface/blocksparsearrayinterface.jl index 0015de0b..d45533cb 100644 --- a/src/blocksparsearrayinterface/blocksparsearrayinterface.jl +++ b/src/blocksparsearrayinterface/blocksparsearrayinterface.jl @@ -146,6 +146,14 @@ end return (I1, to_indices(a, Base.tail(inds), Base.tail(I))...) end +@interface ::AbstractBlockSparseArrayInterface function Base.to_indices( + a, inds, I::Tuple{AbstractArray{Bool},Vararg{Any}} +) + bs1 = to_blockindices(inds[1], I[1]) + I1 = BlockIndices(bs1, blockedunitrange_getindices(inds[1], I[1])) + return (I1, to_indices(a, Base.tail(inds), Base.tail(I))...) +end + # Special case when there is no blocking. @interface ::AbstractBlockSparseArrayInterface function Base.to_indices( a, From afc3f8794d83e466926ef51ba5e4f565b7373080 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 3 Jun 2025 12:47:33 -0400 Subject: [PATCH 03/11] Get logical indexing working --- src/BlockArraysExtensions/BlockArraysExtensions.jl | 7 +++++-- src/abstractblocksparsearray/views.jl | 12 +++++++++--- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/src/BlockArraysExtensions/BlockArraysExtensions.jl b/src/BlockArraysExtensions/BlockArraysExtensions.jl index 3f811347..8bb76189 100644 --- a/src/BlockArraysExtensions/BlockArraysExtensions.jl +++ b/src/BlockArraysExtensions/BlockArraysExtensions.jl @@ -376,8 +376,11 @@ function blockrange(axis::AbstractUnitRange, r::BlockIndexVector) return Block(r):Block(r) end -function blockrange(axis::AbstractUnitRange, r::Vector{<:BlockIndexVector}) - return map(Block, r) +function blockrange( + axis::AbstractUnitRange, + r::BlockVector{<:BlockIndex{1},<:AbstractVector{<:BlockIndexVector}}, +) + return map(Block, blocks(r)) end function blockrange(axis::AbstractUnitRange, r) diff --git a/src/abstractblocksparsearray/views.jl b/src/abstractblocksparsearray/views.jl index f2451988..f7bf776f 100644 --- a/src/abstractblocksparsearray/views.jl +++ b/src/abstractblocksparsearray/views.jl @@ -247,8 +247,9 @@ end # Block slice of the result of slicing `@view a[2:5, 2:5]`. # TODO: Move this to `BlockArraysExtensions`. -const BlockedSlice = BlockSlice{ - <:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}} +const BlockedSlice = Union{ + BlockSlice{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}}}, + BlockIndices{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexVector}}}, } function Base.view( @@ -269,6 +270,11 @@ function BlockArrays.viewblock( ) where {T,N} return viewblock(a, to_tuple(block)...) end + +_block(x) = error("Not implemented.") +_block(x::BlockSlice) = x.block +_block(x::BlockIndices) = x.blocks + # TODO: Define `@interface BlockSparseArrayInterface() viewblock`. function BlockArrays.viewblock( a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{BlockedSlice,N}}}, @@ -279,7 +285,7 @@ function BlockArrays.viewblock( # TODO: Ideally we would use this but it outputs a Vector, # not a range: # return parentindices(a)[dim].block[I[dim]] - return blocks(parentindices(a)[dim].block)[Int(I[dim])] + return blocks(_block(parentindices(a)[dim]))[Int(I[dim])] end return @view parent(a)[brs...] end From fec1e4be9d17f553a69e62e61f668fd2a6d92a68 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 3 Jun 2025 13:09:10 -0400 Subject: [PATCH 04/11] Fix ambiguity error --- src/abstractblocksparsearray/views.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/abstractblocksparsearray/views.jl b/src/abstractblocksparsearray/views.jl index f7bf776f..b3072919 100644 --- a/src/abstractblocksparsearray/views.jl +++ b/src/abstractblocksparsearray/views.jl @@ -271,9 +271,8 @@ function BlockArrays.viewblock( return viewblock(a, to_tuple(block)...) end -_block(x) = error("Not implemented.") -_block(x::BlockSlice) = x.block -_block(x::BlockIndices) = x.blocks +blockedslice_blocks(x::BlockSlice) = x.block +blockedslice_blocks(x::BlockIndices) = x.blocks # TODO: Define `@interface BlockSparseArrayInterface() viewblock`. function BlockArrays.viewblock( @@ -285,7 +284,7 @@ function BlockArrays.viewblock( # TODO: Ideally we would use this but it outputs a Vector, # not a range: # return parentindices(a)[dim].block[I[dim]] - return blocks(_block(parentindices(a)[dim]))[Int(I[dim])] + return blocks(blockedslice_blocks(parentindices(a)[dim]))[Int(I[dim])] end return @view parent(a)[brs...] end From 8a9bc06b3c6d61b528da1e51d1399addca6d1c01 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 3 Jun 2025 13:30:41 -0400 Subject: [PATCH 05/11] Slightly more generally boolean indices --- src/BlockArraysExtensions/blockedunitrange.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/BlockArraysExtensions/blockedunitrange.jl b/src/BlockArraysExtensions/blockedunitrange.jl index eb05e4dc..35884f7f 100644 --- a/src/BlockArraysExtensions/blockedunitrange.jl +++ b/src/BlockArraysExtensions/blockedunitrange.jl @@ -135,7 +135,7 @@ end # TODO: Move this to a `BlockArraysExtensions` library. function blockedunitrange_getindices( - a::AbstractBlockedUnitRange, indices::Vector{<:Integer} + a::AbstractBlockedUnitRange, indices::AbstractVector{<:Integer} ) return map(index -> a[index], indices) end @@ -171,7 +171,7 @@ function blockedunitrange_getindices( end function blockedunitrange_getindices( - a::AbstractBlockedUnitRange, indices::AbstractArray{Bool} + a::AbstractBlockedUnitRange, indices::AbstractVector{Bool} ) blocked_indices = BlockedVector(indices, axes(a)) return mortar( From 5b5d8d4b9a5ac493598196f92802df8eb759dae0 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 3 Jun 2025 17:32:46 -0400 Subject: [PATCH 06/11] Fix for empty blocks --- src/BlockArraysExtensions/blockedunitrange.jl | 24 +++++++++---------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/src/BlockArraysExtensions/blockedunitrange.jl b/src/BlockArraysExtensions/blockedunitrange.jl index 35884f7f..2e916729 100644 --- a/src/BlockArraysExtensions/blockedunitrange.jl +++ b/src/BlockArraysExtensions/blockedunitrange.jl @@ -174,13 +174,12 @@ function blockedunitrange_getindices( a::AbstractBlockedUnitRange, indices::AbstractVector{Bool} ) blocked_indices = BlockedVector(indices, axes(a)) - return mortar( - map(Base.OneTo(blocklength(blocked_indices))) do b - binds = blocked_indices[Block(b)] - bstart = blockfirsts(only(axes(blocked_indices)))[b] - return findall(binds) .+ (bstart - 1) - end, - ) + bs = map(Base.OneTo(blocklength(blocked_indices))) do b + binds = blocked_indices[Block(b)] + bstart = blockfirsts(only(axes(blocked_indices)))[b] + return findall(binds) .+ (bstart - 1) + end + return mortar(filter(!isempty, bs)) end # TODO: Move this to a `BlockArraysExtensions` library. @@ -224,12 +223,11 @@ Base.copy(a::BlockIndexVector) = BlockIndexVector(a.block, copy(a.indices)) function to_blockindices(a::AbstractBlockedUnitRange{<:Integer}, I::AbstractArray{Bool}) I_blocks = blocks(BlockedVector(I, blocklengths(a))) - return mortar( - map(eachindex(I_blocks)) do b - I_b = findall(I_blocks[b]) - BlockIndexVector(Block(b), I_b) - end, - ) + I′_blocks = map(eachindex(I_blocks)) do b + I_b = findall(I_blocks[b]) + BlockIndexVector(Block(b), I_b) + end + return mortar(filter(!isempty, I′_blocks)) end # This handles non-blocked slices. From 199593449a63f4f510c88782a36c9a25142b4a42 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 3 Jun 2025 22:18:48 -0400 Subject: [PATCH 07/11] Fix more slicing operations --- .../BlockArraysExtensions.jl | 5 +-- src/abstractblocksparsearray/views.jl | 34 ++++++++++++++----- 2 files changed, 29 insertions(+), 10 deletions(-) diff --git a/src/BlockArraysExtensions/BlockArraysExtensions.jl b/src/BlockArraysExtensions/BlockArraysExtensions.jl index 8bb76189..1bbfed0b 100644 --- a/src/BlockArraysExtensions/BlockArraysExtensions.jl +++ b/src/BlockArraysExtensions/BlockArraysExtensions.jl @@ -158,8 +158,9 @@ end const BlockSliceCollection = Union{ Base.Slice,BlockSlice{<:BlockRange{1}},BlockIndices{<:Vector{<:Block{1}}} } -const SubBlockSliceCollection = BlockIndices{ - <:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}} +const SubBlockSliceCollection = Union{ + BlockSlice{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}}}, + BlockIndices{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexVector}}}, } # TODO: This is type piracy. This is used in `reindex` when making diff --git a/src/abstractblocksparsearray/views.jl b/src/abstractblocksparsearray/views.jl index b3072919..436666a9 100644 --- a/src/abstractblocksparsearray/views.jl +++ b/src/abstractblocksparsearray/views.jl @@ -92,11 +92,14 @@ end # TODO: Move to `GradedUnitRanges` or `BlockArraysExtensions`. to_block(I::Block{1}) = I to_block(I::BlockIndexRange{1}) = Block(I) +to_block(I::BlockIndexVector) = Block(I) to_block_indices(I::Block{1}) = Colon() to_block_indices(I::BlockIndexRange{1}) = only(I.indices) +to_block_indices(I::BlockIndexVector) = I.indices function Base.view( - a::AbstractBlockSparseArray{<:Any,N}, I::Vararg{Union{Block{1},BlockIndexRange{1}},N} + a::AbstractBlockSparseArray{<:Any,N}, + I::Vararg{Union{Block{1},BlockIndexRange{1},BlockIndexVector},N}, ) where {N} return @views a[to_block.(I)...][to_block_indices.(I)...] end @@ -108,7 +111,7 @@ function Base.view( end function Base.view( a::SubArray{T,N,<:AbstractBlockSparseArray{T,N}}, - I::Vararg{Union{Block{1},BlockIndexRange{1}},N}, + I::Vararg{Union{Block{1},BlockIndexRange{1},BlockIndexVector},N}, ) where {T,N} return @views a[to_block.(I)...][to_block_indices.(I)...] end @@ -205,8 +208,14 @@ function BlockArrays.viewblock( end function to_blockindexrange( - a::BlockIndices{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}}}, - I::Block{1}, + a::BlockSlice{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}}}, I::Block{1} +) + # TODO: Ideally we would just use `a.blocks[I]` but that doesn't + # work right now. + return blocks(a.block)[Int(I)] +end +function to_blockindexrange( + a::BlockIndices{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexVector}}}, I::Block{1} ) # TODO: Ideally we would just use `a.blocks[I]` but that doesn't # work right now. @@ -247,10 +256,13 @@ end # Block slice of the result of slicing `@view a[2:5, 2:5]`. # TODO: Move this to `BlockArraysExtensions`. -const BlockedSlice = Union{ - BlockSlice{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}}}, - BlockIndices{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexVector}}}, +const BlockIndexRangeSlice = BlockSlice{ + <:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}} +} +const BlockIndexVectorSlice = BlockSlice{ + <:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexVector}} } +const BlockedSlice = Union{BlockIndexRangeSlice,BlockIndexVectorSlice} function Base.view( a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{BlockedSlice,N}}}, @@ -258,9 +270,15 @@ function Base.view( ) where {T,N} return viewblock(a, block) end +function Base.view( + a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{BlockIndexRangeSlice,N}}}, + block::Union{Block{N},BlockIndexRange{N}}, +) where {T,N} + return viewblock(a, block) +end function Base.view( a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{BlockedSlice,N}}}, - block::Vararg{Union{Block{1},BlockIndexRange{1}},N}, + block::Vararg{Union{Block{1},BlockIndexRange{1},BlockIndexVector},N}, ) where {T,N} return viewblock(a, block...) end From 77c269c8b3ec49d838732b452727367e9d338b4f Mon Sep 17 00:00:00 2001 From: mtfishman Date: Wed, 4 Jun 2025 10:41:53 -0400 Subject: [PATCH 08/11] Fix tests --- .../BlockArraysExtensions.jl | 12 +++++- src/abstractblocksparsearray/views.jl | 37 +++++++++++-------- src/factorizations/svd.jl | 4 +- test/test_basics.jl | 6 ++- 4 files changed, 40 insertions(+), 19 deletions(-) diff --git a/src/BlockArraysExtensions/BlockArraysExtensions.jl b/src/BlockArraysExtensions/BlockArraysExtensions.jl index 1bbfed0b..5c9f1231 100644 --- a/src/BlockArraysExtensions/BlockArraysExtensions.jl +++ b/src/BlockArraysExtensions/BlockArraysExtensions.jl @@ -158,9 +158,17 @@ end const BlockSliceCollection = Union{ Base.Slice,BlockSlice{<:BlockRange{1}},BlockIndices{<:Vector{<:Block{1}}} } +const BlockIndexRangeSlice = BlockSlice{ + <:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}} +} +const BlockIndexRangeSlices = BlockIndices{ + <:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}} +} +const BlockIndexVectorSlices = BlockIndices{ + <:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexVector}} +} const SubBlockSliceCollection = Union{ - BlockSlice{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}}}, - BlockIndices{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexVector}}}, + BlockIndexRangeSlice,BlockIndexRangeSlices,BlockIndexVectorSlices } # TODO: This is type piracy. This is used in `reindex` when making diff --git a/src/abstractblocksparsearray/views.jl b/src/abstractblocksparsearray/views.jl index 436666a9..a7c93aed 100644 --- a/src/abstractblocksparsearray/views.jl +++ b/src/abstractblocksparsearray/views.jl @@ -214,6 +214,13 @@ function to_blockindexrange( # work right now. return blocks(a.block)[Int(I)] end +function to_blockindexrange( + a::BlockIndices{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange}}}, I::Block{1} +) + # TODO: Ideally we would just use `a.blocks[I]` but that doesn't + # work right now. + return blocks(a.blocks)[Int(I)] +end function to_blockindexrange( a::BlockIndices{<:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexVector}}}, I::Block{1} ) @@ -254,18 +261,10 @@ function BlockArrays.viewblock( return view(viewblock(a, Block.(block)...), map(b -> only(b.indices), block)...) end -# Block slice of the result of slicing `@view a[2:5, 2:5]`. -# TODO: Move this to `BlockArraysExtensions`. -const BlockIndexRangeSlice = BlockSlice{ - <:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexRange{1}}} -} -const BlockIndexVectorSlice = BlockSlice{ - <:BlockVector{<:BlockIndex{1},<:Vector{<:BlockIndexVector}} -} -const BlockedSlice = Union{BlockIndexRangeSlice,BlockIndexVectorSlice} - function Base.view( - a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{BlockedSlice,N}}}, + a::SubArray{ + T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{SubBlockSliceCollection,N}} + }, block::Union{Block{N},BlockIndexRange{N}}, ) where {T,N} return viewblock(a, block) @@ -277,13 +276,17 @@ function Base.view( return viewblock(a, block) end function Base.view( - a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{BlockedSlice,N}}}, + a::SubArray{ + T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{SubBlockSliceCollection,N}} + }, block::Vararg{Union{Block{1},BlockIndexRange{1},BlockIndexVector},N}, ) where {T,N} return viewblock(a, block...) end function BlockArrays.viewblock( - a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{BlockedSlice,N}}}, + a::SubArray{ + T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{SubBlockSliceCollection,N}} + }, block::Union{Block{N},BlockIndexRange{N}}, ) where {T,N} return viewblock(a, to_tuple(block)...) @@ -294,7 +297,9 @@ blockedslice_blocks(x::BlockIndices) = x.blocks # TODO: Define `@interface BlockSparseArrayInterface() viewblock`. function BlockArrays.viewblock( - a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{BlockedSlice,N}}}, + a::SubArray{ + T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{SubBlockSliceCollection,N}} + }, I::Vararg{Block{1},N}, ) where {T,N} # TODO: Use `reindex`, `to_indices`, etc. @@ -308,7 +313,9 @@ function BlockArrays.viewblock( end # TODO: Define `@interface BlockSparseArrayInterface() viewblock`. function BlockArrays.viewblock( - a::SubArray{T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{BlockedSlice,N}}}, + a::SubArray{ + T,N,<:AbstractBlockSparseArray{T,N},<:Tuple{Vararg{SubBlockSliceCollection,N}} + }, block::Vararg{BlockIndexRange{1},N}, ) where {T,N} return view(viewblock(a, Block.(block)...), map(b -> only(b.indices), block)...) diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl index c73b8e24..635da79f 100644 --- a/src/factorizations/svd.jl +++ b/src/factorizations/svd.jl @@ -12,7 +12,9 @@ struct BlockPermutedDiagonalAlgorithm{A<:MatrixAlgebraKit.AbstractAlgorithm} <: alg::A end -function MatrixAlgebraKit.default_svd_algorithm(A::AbstractBlockSparseMatrix; kwargs...) +function MatrixAlgebraKit.default_svd_algorithm( + A::Type{<:AbstractBlockSparseMatrix}; kwargs... +) blocktype(A) <: StridedMatrix{<:LinearAlgebra.BLAS.BlasFloat} || error("unsupported type: $(blocktype(A))") # TODO: this is a hardcoded for now to get around this function not being defined in the diff --git a/test/test_basics.jl b/test/test_basics.jl index b0ee9a58..d8321aba 100644 --- a/test/test_basics.jl +++ b/test/test_basics.jl @@ -588,7 +588,11 @@ arrayts = (Array, JLArray) @views for b in [Block(1, 2), Block(2, 1)] a[b] = dev(randn(elt, size(a[b]))) end - for b in (a[2:4, 2:4], @view(a[2:4, 2:4])) + I = 2:4 + I′, J′ = falses.(size(a)) + I′[I] .= true + J′[I] .= true + for b in (a[I, I], @view(a[I, I]), a[I′, J′], @view(a[I′, J′])) @allowscalar @test b == Array(a)[2:4, 2:4] @test size(b) == (3, 3) @test blocksize(b) == (2, 2) From f0e3ea2ee3a53d476a23bd730f86afc7e3516d63 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Wed, 4 Jun 2025 12:30:06 -0400 Subject: [PATCH 09/11] Fix factorizations --- src/factorizations/lq.jl | 19 +++++++++++-------- src/factorizations/polar.jl | 18 ++++++++++++++---- src/factorizations/qr.jl | 11 +++++++---- src/factorizations/svd.jl | 7 +++++-- 4 files changed, 37 insertions(+), 18 deletions(-) diff --git a/src/factorizations/lq.jl b/src/factorizations/lq.jl index 4a07cfa6..980b7f70 100644 --- a/src/factorizations/lq.jl +++ b/src/factorizations/lq.jl @@ -1,22 +1,25 @@ -using MatrixAlgebraKit: MatrixAlgebraKit, lq_compact!, lq_full! +using MatrixAlgebraKit: MatrixAlgebraKit, default_lq_algorithm, lq_compact!, lq_full! -# TODO: this is a hardcoded for now to get around this function not being defined in the -# type domain -function default_blocksparse_lq_algorithm(A::AbstractMatrix; kwargs...) +function MatrixAlgebraKit.default_lq_algorithm(A::AbstractBlockSparseMatrix; kwargs...) + return default_lq_algorithm(typeof(A); kwargs...) +end +function MatrixAlgebraKit.default_lq_algorithm( + A::Type{<:AbstractBlockSparseMatrix}; kwargs... +) blocktype(A) <: StridedMatrix{<:LinearAlgebra.BLAS.BlasFloat} || error("unsupported type: $(blocktype(A))") alg = MatrixAlgebraKit.LAPACK_HouseholderLQ(; kwargs...) return BlockPermutedDiagonalAlgorithm(alg) end function MatrixAlgebraKit.default_algorithm( - ::typeof(lq_compact!), A::AbstractBlockSparseMatrix; kwargs... + ::typeof(lq_compact!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) - return default_blocksparse_lq_algorithm(A; kwargs...) + return default_lq_algorithm(A; kwargs...) end function MatrixAlgebraKit.default_algorithm( - ::typeof(lq_full!), A::AbstractBlockSparseMatrix; kwargs... + ::typeof(lq_full!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) - return default_blocksparse_lq_algorithm(A; kwargs...) + return default_q_algorithm(A; kwargs...) end function similar_output( diff --git a/src/factorizations/polar.jl b/src/factorizations/polar.jl index 9b9c2831..2a6313b4 100644 --- a/src/factorizations/polar.jl +++ b/src/factorizations/polar.jl @@ -48,12 +48,22 @@ function MatrixAlgebraKit.right_polar!(A::AbstractBlockSparseMatrix, alg::PolarV end function MatrixAlgebraKit.default_algorithm( - ::typeof(left_polar!), a::AbstractBlockSparseMatrix; kwargs... + ::typeof(left_polar!), A::AbstractBlockSparseMatrix; kwargs... ) - return PolarViaSVD(default_algorithm(svd_compact!, a; kwargs...)) + return default_algorithm(left_polar!, typeof(A); kwargs...) end function MatrixAlgebraKit.default_algorithm( - ::typeof(right_polar!), a::AbstractBlockSparseMatrix; kwargs... + ::typeof(left_polar!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) - return PolarViaSVD(default_algorithm(svd_compact!, a; kwargs...)) + return PolarViaSVD(default_algorithm(svd_compact!, A; kwargs...)) +end +function MatrixAlgebraKit.default_algorithm( + ::typeof(right_polar!), A::AbstractBlockSparseMatrix; kwargs... +) + return default_algorithm(right_polar!, typeof(A); kwargs...) +end +function MatrixAlgebraKit.default_algorithm( + ::typeof(right_polar!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... +) + return PolarViaSVD(default_algorithm(svd_compact!, A; kwargs...)) end diff --git a/src/factorizations/qr.jl b/src/factorizations/qr.jl index 55a0b93e..775b7094 100644 --- a/src/factorizations/qr.jl +++ b/src/factorizations/qr.jl @@ -1,21 +1,24 @@ using MatrixAlgebraKit: MatrixAlgebraKit, default_qr_algorithm, lq_compact!, lq_full!, qr_compact!, qr_full! -# TODO: this is a hardcoded for now to get around this function not being defined in the -# type domain function MatrixAlgebraKit.default_qr_algorithm(A::AbstractBlockSparseMatrix; kwargs...) + return default_qr_algorithm(typeof(A); kwargs...) +end +function MatrixAlgebraKit.default_qr_algorithm( + A::Type{<:AbstractBlockSparseMatrix}; kwargs... +) blocktype(A) <: StridedMatrix{<:LinearAlgebra.BLAS.BlasFloat} || error("unsupported type: $(blocktype(A))") alg = MatrixAlgebraKit.LAPACK_HouseholderQR(; kwargs...) return BlockPermutedDiagonalAlgorithm(alg) end function MatrixAlgebraKit.default_algorithm( - ::typeof(qr_compact!), A::AbstractBlockSparseMatrix; kwargs... + ::typeof(qr_compact!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) return default_qr_algorithm(A; kwargs...) end function MatrixAlgebraKit.default_algorithm( - ::typeof(qr_full!), A::AbstractBlockSparseMatrix; kwargs... + ::typeof(qr_full!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) return default_qr_algorithm(A; kwargs...) end diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl index 635da79f..a3ab06ee 100644 --- a/src/factorizations/svd.jl +++ b/src/factorizations/svd.jl @@ -12,6 +12,9 @@ struct BlockPermutedDiagonalAlgorithm{A<:MatrixAlgebraKit.AbstractAlgorithm} <: alg::A end +function MatrixAlgebraKit.default_svd_algorithm(A::AbstractBlockSparseMatrix; kwargs...) + return default_svd_algorithm(typeof(A), kwargs...) +end function MatrixAlgebraKit.default_svd_algorithm( A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) @@ -25,12 +28,12 @@ function MatrixAlgebraKit.default_svd_algorithm( end function MatrixAlgebraKit.default_algorithm( - f::typeof(svd_compact!), A::AbstractBlockSparseMatrix; kwargs... + f::typeof(svd_compact!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) return default_svd_algorithm(A; kwargs...) end function MatrixAlgebraKit.default_algorithm( - f::typeof(svd_full!), A::AbstractBlockSparseMatrix; kwargs... + f::typeof(svd_full!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) return default_svd_algorithm(A; kwargs...) end From a5e0c3f298ad74ca7ccfbb16ac0b40173c0018e4 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Wed, 4 Jun 2025 12:53:52 -0400 Subject: [PATCH 10/11] Try fixing tests --- .../wrappedabstractblocksparsearray.jl | 7 +++ src/factorizations/lq.jl | 22 +++++---- src/factorizations/polar.jl | 49 +++++++++++-------- src/factorizations/qr.jl | 20 +++++--- src/factorizations/svd.jl | 31 +++++++----- 5 files changed, 80 insertions(+), 49 deletions(-) diff --git a/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl b/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl index 91705772..c22907b8 100644 --- a/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl +++ b/src/abstractblocksparsearray/wrappedabstractblocksparsearray.jl @@ -45,6 +45,13 @@ function Base.to_indices( ) return @interface BlockSparseArrayInterface() to_indices(a, inds, I) end +# Fix ambiguity error with Base for logical indexing in Julia 1.10. +# TODO: Delete this once we drop support for Julia 1.10. +function Base.to_indices( + a::AnyAbstractBlockSparseArray, inds, I::Union{Tuple{BitArray{N}},Tuple{Array{Bool,N}}} +) where {N} + return @interface BlockSparseArrayInterface() to_indices(a, inds, I) +end # a[[Block(2), Block(1)], [Block(2), Block(1)]] function Base.to_indices( diff --git a/src/factorizations/lq.jl b/src/factorizations/lq.jl index 980b7f70..acb1046f 100644 --- a/src/factorizations/lq.jl +++ b/src/factorizations/lq.jl @@ -1,25 +1,29 @@ using MatrixAlgebraKit: MatrixAlgebraKit, default_lq_algorithm, lq_compact!, lq_full! +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. function MatrixAlgebraKit.default_lq_algorithm(A::AbstractBlockSparseMatrix; kwargs...) return default_lq_algorithm(typeof(A); kwargs...) end -function MatrixAlgebraKit.default_lq_algorithm( - A::Type{<:AbstractBlockSparseMatrix}; kwargs... -) - blocktype(A) <: StridedMatrix{<:LinearAlgebra.BLAS.BlasFloat} || - error("unsupported type: $(blocktype(A))") - alg = MatrixAlgebraKit.LAPACK_HouseholderLQ(; kwargs...) - return BlockPermutedDiagonalAlgorithm(alg) -end + +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. function MatrixAlgebraKit.default_algorithm( ::typeof(lq_compact!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) return default_lq_algorithm(A; kwargs...) end + +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. function MatrixAlgebraKit.default_algorithm( ::typeof(lq_full!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) - return default_q_algorithm(A; kwargs...) + return default_lq_algorithm(A; kwargs...) +end + +function MatrixAlgebraKit.default_lq_algorithm( + A::Type{<:AbstractBlockSparseMatrix}; kwargs... +) + alg = default_lq_algorithm(blocktype(A); kwargs...) + return BlockPermutedDiagonalAlgorithm(alg) end function similar_output( diff --git a/src/factorizations/polar.jl b/src/factorizations/polar.jl index 2a6313b4..d8abba7c 100644 --- a/src/factorizations/polar.jl +++ b/src/factorizations/polar.jl @@ -7,6 +7,34 @@ using MatrixAlgebraKit: right_polar!, svd_compact! +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. +function MatrixAlgebraKit.default_algorithm( + ::typeof(left_polar!), A::AbstractBlockSparseMatrix; kwargs... +) + return default_algorithm(left_polar!, typeof(A); kwargs...) +end + +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. +function MatrixAlgebraKit.default_algorithm( + ::typeof(left_polar!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... +) + return PolarViaSVD(default_algorithm(svd_compact!, A; kwargs...)) +end + +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. +function MatrixAlgebraKit.default_algorithm( + ::typeof(right_polar!), A::AbstractBlockSparseMatrix; kwargs... +) + return default_algorithm(right_polar!, typeof(A); kwargs...) +end + +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. +function MatrixAlgebraKit.default_algorithm( + ::typeof(right_polar!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... +) + return PolarViaSVD(default_algorithm(svd_compact!, A; kwargs...)) +end + function MatrixAlgebraKit.check_input(::typeof(left_polar!), A::AbstractBlockSparseMatrix) @views for I in eachblockstoredindex(A) m, n = size(A[I]) @@ -46,24 +74,3 @@ function MatrixAlgebraKit.right_polar!(A::AbstractBlockSparseMatrix, alg::PolarV P = U * S * copy(U') return (P, Wᴴ) end - -function MatrixAlgebraKit.default_algorithm( - ::typeof(left_polar!), A::AbstractBlockSparseMatrix; kwargs... -) - return default_algorithm(left_polar!, typeof(A); kwargs...) -end -function MatrixAlgebraKit.default_algorithm( - ::typeof(left_polar!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... -) - return PolarViaSVD(default_algorithm(svd_compact!, A; kwargs...)) -end -function MatrixAlgebraKit.default_algorithm( - ::typeof(right_polar!), A::AbstractBlockSparseMatrix; kwargs... -) - return default_algorithm(right_polar!, typeof(A); kwargs...) -end -function MatrixAlgebraKit.default_algorithm( - ::typeof(right_polar!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... -) - return PolarViaSVD(default_algorithm(svd_compact!, A; kwargs...)) -end diff --git a/src/factorizations/qr.jl b/src/factorizations/qr.jl index 775b7094..74e587b7 100644 --- a/src/factorizations/qr.jl +++ b/src/factorizations/qr.jl @@ -1,28 +1,32 @@ using MatrixAlgebraKit: MatrixAlgebraKit, default_qr_algorithm, lq_compact!, lq_full!, qr_compact!, qr_full! +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. function MatrixAlgebraKit.default_qr_algorithm(A::AbstractBlockSparseMatrix; kwargs...) return default_qr_algorithm(typeof(A); kwargs...) end -function MatrixAlgebraKit.default_qr_algorithm( - A::Type{<:AbstractBlockSparseMatrix}; kwargs... -) - blocktype(A) <: StridedMatrix{<:LinearAlgebra.BLAS.BlasFloat} || - error("unsupported type: $(blocktype(A))") - alg = MatrixAlgebraKit.LAPACK_HouseholderQR(; kwargs...) - return BlockPermutedDiagonalAlgorithm(alg) -end + +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. function MatrixAlgebraKit.default_algorithm( ::typeof(qr_compact!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) return default_qr_algorithm(A; kwargs...) end + +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. function MatrixAlgebraKit.default_algorithm( ::typeof(qr_full!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) return default_qr_algorithm(A; kwargs...) end +function MatrixAlgebraKit.default_qr_algorithm( + A::Type{<:AbstractBlockSparseMatrix}; kwargs... +) + alg = default_qr_algorithm(blocktype(A); kwargs...) + return BlockPermutedDiagonalAlgorithm(alg) +end + function similar_output( ::typeof(qr_compact!), A, R_axis, alg::MatrixAlgebraKit.AbstractAlgorithm ) diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl index a3ab06ee..1d4d88bf 100644 --- a/src/factorizations/svd.jl +++ b/src/factorizations/svd.jl @@ -1,5 +1,14 @@ using MatrixAlgebraKit: MatrixAlgebraKit, default_svd_algorithm, svd_compact!, svd_full! +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. +using MatrixAlgebraKit: TruncatedAlgorithm, select_truncation, svd_trunc! +function MatrixAlgebraKit.select_algorithm( + ::typeof(svd_trunc!), A::Type{<:AbstractBlockSparseMatrix}, alg; trunc=nothing, kwargs... +) + alg_svd = select_algorithm(svd_compact!, A, alg; kwargs...) + return TruncatedAlgorithm(alg_svd, select_truncation(trunc)) +end + """ BlockPermutedDiagonalAlgorithm(A::MatrixAlgebraKit.AbstractAlgorithm) @@ -12,32 +21,32 @@ struct BlockPermutedDiagonalAlgorithm{A<:MatrixAlgebraKit.AbstractAlgorithm} <: alg::A end +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. function MatrixAlgebraKit.default_svd_algorithm(A::AbstractBlockSparseMatrix; kwargs...) return default_svd_algorithm(typeof(A), kwargs...) end -function MatrixAlgebraKit.default_svd_algorithm( - A::Type{<:AbstractBlockSparseMatrix}; kwargs... -) - blocktype(A) <: StridedMatrix{<:LinearAlgebra.BLAS.BlasFloat} || - error("unsupported type: $(blocktype(A))") - # TODO: this is a hardcoded for now to get around this function not being defined in the - # type domain - # alg = MatrixAlgebraKit.default_algorithm(f, blocktype(A); kwargs...) - alg = MatrixAlgebraKit.LAPACK_DivideAndConquer(; kwargs...) - return BlockPermutedDiagonalAlgorithm(alg) -end +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. function MatrixAlgebraKit.default_algorithm( f::typeof(svd_compact!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) return default_svd_algorithm(A; kwargs...) end + +# TODO: Delete once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/32 is merged. function MatrixAlgebraKit.default_algorithm( f::typeof(svd_full!), A::Type{<:AbstractBlockSparseMatrix}; kwargs... ) return default_svd_algorithm(A; kwargs...) end +function MatrixAlgebraKit.default_svd_algorithm( + A::Type{<:AbstractBlockSparseMatrix}; kwargs... +) + alg = default_svd_algorithm(blocktype(A); kwargs...) + return BlockPermutedDiagonalAlgorithm(alg) +end + function similar_output( ::typeof(svd_compact!), A, S_axes, alg::MatrixAlgebraKit.AbstractAlgorithm ) From 2095a54ee70f00f2ea0ed5430ee0edc9004cc522 Mon Sep 17 00:00:00 2001 From: Matt Fishman Date: Wed, 4 Jun 2025 13:17:36 -0400 Subject: [PATCH 11/11] Bump version properly --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ef35add8..327cea0e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BlockSparseArrays" uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" authors = ["ITensor developers and contributors"] -version = "0.6.10" +version = "0.6.9" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"