More generalizations for generic block types (#144)

Author: mtfishman

Project.toml

@@ -1,7 +1,7 @@
 name = "BlockSparseArrays"
 uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
 authors = ["ITensor developers <[email protected]> and contributors"]
-version = "0.7.9"
+version = "0.7.10"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -44,7 +44,7 @@ SparseArraysBase = "0.5"
 SplitApplyCombine = "1.2.3"
 TensorAlgebra = "0.3.2"
 Test = "1.10"
-TypeParameterAccessors = "0.4"
+TypeParameterAccessors = "0.4.1"
 julia = "1.10"
 
 [extras]

src/blocksparsearray/blocksparsearray.jl

@@ -9,6 +9,7 @@ using BlockArrays:
 using DerivableInterfaces: @interface
 using Dictionaries: Dictionary
 using SparseArraysBase: SparseArrayDOK
+using TypeParameterAccessors: similartype
 
 """
     SparseArrayDOK{T}(undef_blocks, axes)
@@ -173,7 +174,9 @@ end
 function BlockSparseArray{T,N}(
   ::UndefInitializer, axes::Tuple{Vararg{AbstractUnitRange{<:Integer},N}}
 ) where {T,N}
-  return BlockSparseArray{T,N,Array{T,N}}(undef, axes)
+  axt = Tuple{blockaxistype.(axes)...}
+  A = similartype(Array{T}, axt)
+  return BlockSparseArray{T,N,A}(undef, axes)
 end
 
 function BlockSparseArray{T,N}(

src/blocksparsearrayinterface/arraylayouts.jl

@@ -11,15 +11,18 @@ using LinearAlgebra: LinearAlgebra, dot, mul!
   return a_dest
 end
 
+function DerivableInterfaces.interface(m::MulAdd)
+  return interface(m.A, m.B, m.C)
+end
+
 function ArrayLayouts.materialize!(
   m::MatMulMatAdd{
     <:BlockLayout{<:SparseLayout},
     <:BlockLayout{<:SparseLayout},
     <:BlockLayout{<:SparseLayout},
   },
 )
-  α, a1, a2, β, a_dest = m.α, m.A, m.B, m.β, m.C
-  @interface BlockSparseArrayInterface() muladd!(m.α, m.A, m.B, m.β, m.C)
+  @interface interface(m) muladd!(m.α, m.A, m.B, m.β, m.C)
   return m.C
 end
 function ArrayLayouts.materialize!(
@@ -29,7 +32,7 @@ function ArrayLayouts.materialize!(
     <:BlockLayout{<:SparseLayout},
   },
 )
-  @interface BlockSparseArrayInterface() matmul!(m)
+  @interface interface(m) matmul!(m)
   return m.C
 end
 
@@ -42,5 +45,5 @@ end
 end
 
 function Base.copy(d::Dot{<:BlockLayout{<:SparseLayout},<:BlockLayout{<:SparseLayout}})
-  return @interface BlockSparseArrayInterface() dot(d.A, d.B)
+  return @interface interface(d.A, d.B) dot(d.A, d.B)
 end

src/factorizations/svd.jl

@@ -1,5 +1,7 @@
+using DiagonalArrays: diagonaltype
 using MatrixAlgebraKit:
   MatrixAlgebraKit, check_input, default_svd_algorithm, svd_compact!, svd_full!
+using TypeParameterAccessors: realtype
 
 """
     BlockPermutedDiagonalAlgorithm(A::MatrixAlgebraKit.AbstractAlgorithm)
@@ -24,10 +26,7 @@ function similar_output(
   ::typeof(svd_compact!), A, S_axes, alg::MatrixAlgebraKit.AbstractAlgorithm
 )
   U = similar(A, axes(A, 1), S_axes[1])
-  T = real(eltype(A))
-  # TODO: this should be replaced with a more general similar function that can handle setting
-  # the blocktype and element type - something like S = similar(A, BlockType(...))
-  S = BlockSparseMatrix{T,Diagonal{T,Vector{T}}}(undef, S_axes)
+  S = similar(A, BlockType(diagonaltype(realtype(blocktype(A)))), S_axes)
   Vt = similar(A, S_axes[2], axes(A, 2))
   return U, S, Vt
 end