Add MatrixAlgebraKit factorizations

Author: lkdvos

Project.toml

@@ -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"

src/BlockSparseArrays.jl

@@ -8,8 +8,6 @@ export BlockSparseArray,
   eachstoredblock,
   sparsemortar
 
-# factorizations
-include("factorizations/svd.jl")
 
 # possible upstream contributions
 include("BlockArraysExtensions/blockedunitrange.jl")
@@ -45,4 +43,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 = similar(V, b_size[1])
+    return LinearAlgebra.Diagonal{T,V}(diag)
+  else
+    return zeros(T, b_size...)
+  end
+end