Add MatrixAlgebraKit svd

Project.toml

@@ -1,7 +1,7 @@
 name = "BlockSparseArrays"
 uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
 authors = ["ITensor developers <[email protected]> and contributors"]
-version = "0.4.6"
+version = "0.5.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -15,6 +15,7 @@ GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 MapBroadcast = "ebd9b9da-f48d-417c-9660-449667d60261"
+MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
 SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208"
 SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
@@ -38,6 +39,7 @@ GPUArraysCore = "0.1.0, 0.2"
 LinearAlgebra = "1.10"
 MacroTools = "0.5.13"
 MapBroadcast = "0.1.5"
+MatrixAlgebraKit = "0.1.2"
 SparseArraysBase = "0.5"
 SplitApplyCombine = "1.2.3"
 TensorAlgebra = "0.3.2"

docs/Project.toml

@@ -6,6 +6,6 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 
 [compat]
 BlockArrays = "1"
-BlockSparseArrays = "0.4"
+BlockSparseArrays = "0.5"
 Documenter = "1"
 Literate = "2"

examples/Project.toml

@@ -5,5 +5,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 BlockArrays = "1"
-BlockSparseArrays = "0.4"
+BlockSparseArrays = "0.5"
 Test = "1"

src/BlockSparseArrays.jl

@@ -8,9 +8,6 @@ export BlockSparseArray,
   eachstoredblock,
   sparsemortar
 
-# factorizations
-include("factorizations/svd.jl")
-
 # possible upstream contributions
 include("BlockArraysExtensions/blockedunitrange.jl")
 include("BlockArraysExtensions/BlockArraysExtensions.jl")
@@ -45,4 +42,7 @@ include("blocksparsearray/blockdiagonalarray.jl")
 
 include("BlockArraysSparseArraysBaseExt/BlockArraysSparseArraysBaseExt.jl")
 
+# factorizations
+include("factorizations/svd.jl")
+
 end

src/blocksparsearrayinterface/getunstoredblock.jl

@@ -23,3 +23,18 @@ end
 ) where {N}
   return f(a, Tuple(I)...)
 end
+
+# TODO: this is a hack and is also type-unstable
+function (f::GetUnstoredBlock)(
+  a::AbstractMatrix{LinearAlgebra.Diagonal{T,V}}, I::Vararg{Int,2}
+) where {T,V}
+  b_size = ntuple(ndims(a)) do d
+    return length(f.axes[d][Block(I[d])])
+  end
+  if I[1] == I[2]
+    diag = zero!(similar(V, b_size[1]))
+    return LinearAlgebra.Diagonal{T,V}(diag)
+  else
+    return zeros(T, b_size...)
+  end
+end