diff --git a/Project.toml b/Project.toml index 327cea0e..b0a350bc 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.7.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/docs/Project.toml b/docs/Project.toml index d357726b..da8df1ab 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -6,6 +6,6 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" [compat] BlockArrays = "1" -BlockSparseArrays = "0.6" +BlockSparseArrays = "0.7" Documenter = "1" Literate = "2" diff --git a/examples/Project.toml b/examples/Project.toml index d5e5b642..b7d93053 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -5,5 +5,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] BlockArrays = "1" -BlockSparseArrays = "0.6" +BlockSparseArrays = "0.7" Test = "1" diff --git a/src/blocksparsearrayinterface/blocksparsearrayinterface.jl b/src/blocksparsearrayinterface/blocksparsearrayinterface.jl index d45533cb..d11f604c 100644 --- a/src/blocksparsearrayinterface/blocksparsearrayinterface.jl +++ b/src/blocksparsearrayinterface/blocksparsearrayinterface.jl @@ -432,9 +432,11 @@ function Base.isassigned(a::SparseSubArrayBlocks{<:Any,N}, I::Vararg{Int,N}) whe end function SparseArraysBase.eachstoredindex(::IndexCartesian, a::SparseSubArrayBlocks) - return filter(eachindex(a)) do I + isempty(a) && return CartesianIndex{ndims(a)}[] + inds = filter(eachindex(a)) do I return isstored(a, I) end + return inds ## # TODO: This only works for blockwise slices, i.e. slices using ## # `BlockSliceCollection`. diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index e3362128..175e78f3 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -45,69 +45,11 @@ function MatrixAlgebraKit.findtruncated( return indexmask end -function similar_truncate( - ::typeof(svd_trunc!), - (U, S, Vᴴ)::TBlockUSVᴴ, - strategy::BlockPermutedDiagonalTruncationStrategy, - indexmask=MatrixAlgebraKit.findtruncated(diagview(S), strategy), -) - ax = axes(S, 1) - counter = Base.Fix1(count, Base.Fix1(getindex, indexmask)) - s_lengths = filter!(>(0), map(counter, blocks(ax))) - s_axis = blockedrange(s_lengths) - Ũ = similar(U, axes(U, 1), s_axis) - S̃ = similar(S, s_axis, s_axis) - Ṽᴴ = similar(Vᴴ, s_axis, axes(Vᴴ, 2)) - return Ũ, S̃, Ṽᴴ -end - function MatrixAlgebraKit.truncate!( ::typeof(svd_trunc!), (U, S, Vᴴ)::TBlockUSVᴴ, strategy::BlockPermutedDiagonalTruncationStrategy, ) - indexmask = MatrixAlgebraKit.findtruncated(diagview(S), strategy) - - # first determine the block structure of the output to avoid having assumptions on the - # data structures - Ũ, S̃, Ṽᴴ = similar_truncate(svd_trunc!, (U, S, Vᴴ), strategy, indexmask) - - # then loop over the blocks and assign the data - # TODO: figure out if we can presort and loop over the blocks - - # for now this has issues with missing blocks - bI_Us = collect(eachblockstoredindex(U)) - bI_Ss = collect(eachblockstoredindex(S)) - bI_Vᴴs = collect(eachblockstoredindex(Vᴴ)) - - I′ = 0 # number of skipped blocks that got fully truncated - ax = axes(S, 1) - for I in 1:blocksize(ax, 1) - b = ax[Block(I)] - mask = indexmask[b] - - if !any(mask) - I′ += 1 - continue - end - - bU_id = @something findfirst(x -> last(Tuple(x)) == Block(I), bI_Us) error( - "No U-block found for $I" - ) - bU = Tuple(bI_Us[bU_id]) - Ũ[bU[1], bU[2] - Block(I′)] = view(U, bU...)[:, mask] - - bVᴴ_id = @something findfirst(x -> first(Tuple(x)) == Block(I), bI_Vᴴs) error( - "No Vᴴ-block found for $I" - ) - bVᴴ = Tuple(bI_Vᴴs[bVᴴ_id]) - Ṽᴴ[bVᴴ[1] - Block(I′), bVᴴ[2]] = view(Vᴴ, bVᴴ...)[mask, :] - - bS_id = findfirst(x -> last(Tuple(x)) == Block(I), bI_Ss) - if !isnothing(bS_id) - bS = Tuple(bI_Ss[bS_id]) - S̃[(bS .- Block(I′))...] = Diagonal(diagview(view(S, bS...))[mask]) - end - end - - return Ũ, S̃, Ṽᴴ + I = MatrixAlgebraKit.findtruncated(diagview(S), strategy) + return (U[:, I], S[I, I], Vᴴ[I, :]) end diff --git a/test/Project.toml b/test/Project.toml index e40faef7..66cf5f46 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -23,7 +23,7 @@ Adapt = "4" Aqua = "0.8" ArrayLayouts = "1" BlockArrays = "1" -BlockSparseArrays = "0.6" +BlockSparseArrays = "0.7" DiagonalArrays = "0.3" GPUArraysCore = "0.2" JLArrays = "0.2"