Fix truncated SVD (#39)

Author: mtfishman

Project.toml

@@ -1,7 +1,7 @@
 name = "GradedArrays"
 uuid = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2"
 authors = ["ITensor developers <[email protected]> and contributors"]
-version = "0.4.6"
+version = "0.4.7"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
@@ -25,7 +25,7 @@ GradedArraysTensorAlgebraExt = "TensorAlgebra"
 
 [compat]
 BlockArrays = "1.6.0"
-BlockSparseArrays = "0.6.5"
+BlockSparseArrays = "0.7.0"
 Compat = "4.16.0"
 DerivableInterfaces = "0.4.4"
 FillArrays = "1.13.0"

src/factorizations.jl

@@ -42,29 +42,6 @@ end
 
 const TGradedUSVᴴ = Tuple{<:GradedMatrix,<:GradedMatrix,<:GradedMatrix}
 
-function BlockSparseArrays.similar_truncate(
-  ::typeof(svd_trunc!),
-  (U, S, Vᴴ)::TGradedUSVᴴ,
-  strategy::BlockPermutedDiagonalTruncationStrategy,
-  indexmask=MatrixAlgebraKit.findtruncated(diagview(S), strategy),
-)
-  u_axis, v_axis = axes(S)
-  counter = Base.Fix1(count, Base.Fix1(getindex, indexmask))
-  s_lengths = map(counter, blocks(u_axis))
-  u_sectors = sectors(u_axis) .=> s_lengths
-  v_sectors = sectors(v_axis) .=> s_lengths
-  u_sectors_filtered = filter(>(0) ∘ last, u_sectors)
-  v_sectors_filtered = filter(>(0) ∘ last, v_sectors)
-  u_axis′ = gradedrange(u_sectors_filtered)
-  u_axis = isdual(u_axis) ? dual(u_axis′) : u_axis′
-  v_axis′ = gradedrange(v_sectors_filtered)
-  v_axis = isdual(v_axis) ? dual(v_axis′) : v_axis′
-  Ũ = similar(U, axes(U, 1), dual(u_axis))
-  S̃ = similar(S, u_axis, v_axis)
-  Ṽᴴ = similar(Vᴴ, dual(v_axis), axes(Vᴴ, 2))
-  return Ũ, S̃, Ṽᴴ
-end
-
 function BlockSparseArrays.similar_output(
   ::typeof(qr_compact!), A::GradedMatrix, R_axis, alg::BlockPermutedDiagonalAlgorithm
 )

src/gradedunitrange.jl

@@ -320,3 +320,21 @@ function BlockSparseArrays.blockedunitrange_getindices(
   new_range = blockedrange(new_first, length.(new_axes))
   return GradedUnitRange(new_axes, new_range)
 end
+
+using BlockArrays: BlockedVector
+function BlockSparseArrays.blockedunitrange_getindices(
+  a::AbstractGradedUnitRange, indices::AbstractVector{Bool}
+)
+  blocked_indices = BlockedVector(indices, axes(a))
+  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
+  keep = map(!isempty, bs)
+  secs = sectors(a)[keep]
+  bs = bs[keep]
+  r = gradedrange(secs .=> length.(bs); isdual=isdual(a))
+  I = mortar(bs, (r,))
+  return I
+end

test/Project.toml

@@ -17,7 +17,7 @@ TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 [compat]
 Aqua = "0.8.11"
 BlockArrays = "1.6.0"
-BlockSparseArrays = "0.6"
+BlockSparseArrays = "0.7.0"
 GradedArrays = "0.4"
 LinearAlgebra = "1.10.0"
 MatrixAlgebraKit = "0.2"