Catch more cases of sparse matmul involving adjoints (#60)

Author: mtfishman

Project.toml

@@ -1,7 +1,7 @@
 name = "SparseArraysBase"
 uuid = "0d5efcca-f356-4864-8770-e1ed8d78f208"
 authors = ["ITensor developers <[email protected]> and contributors"]
-version = "0.5.8"
+version = "0.5.9"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

src/abstractsparsearray.jl

@@ -152,3 +152,43 @@ function sparserand!(
     A[I] = v
   end
 end
+
+# Catch some cases that aren't getting caught by the current
+# DerivableInterfaces.jl logic.
+# TODO: Make this more systematic once DerivableInterfaces.jl
+# is rewritten.
+using ArrayLayouts: ArrayLayouts, MemoryLayout
+using LinearAlgebra: LinearAlgebra, Adjoint
+function ArrayLayouts.MemoryLayout(::Type{Transpose{T,P}}) where {T,P<:AbstractSparseMatrix}
+  return MemoryLayout(P)
+end
+function ArrayLayouts.MemoryLayout(::Type{Adjoint{T,P}}) where {T,P<:AbstractSparseMatrix}
+  return MemoryLayout(P)
+end
+function LinearAlgebra.mul!(
+  dest::AbstractMatrix,
+  A::Adjoint{<:Any,<:AbstractSparseMatrix},
+  B::AbstractSparseMatrix,
+  α::Number,
+  β::Number,
+)
+  return ArrayLayouts.mul!(dest, A, B, α, β)
+end
+function LinearAlgebra.mul!(
+  dest::AbstractMatrix,
+  A::AbstractSparseMatrix,
+  B::Adjoint{<:Any,<:AbstractSparseMatrix},
+  α::Number,
+  β::Number,
+)
+  return ArrayLayouts.mul!(dest, A, B, α, β)
+end
+function LinearAlgebra.mul!(
+  dest::AbstractMatrix,
+  A::Adjoint{<:Any,<:AbstractSparseMatrix},
+  B::Adjoint{<:Any,<:AbstractSparseMatrix},
+  α::Number,
+  β::Number,
+)
+  return ArrayLayouts.mul!(dest, A, B, α, β)
+end

test/test_linalg.jl

@@ -16,13 +16,22 @@ const rng = StableRNG(123)
     A = sparserand(rng, T, szA; density)
     B = sparserand(rng, T, szB; density)
 
-    check1 = mul!(Array(C), Array(A), Array(B))
-    @test mul!(copy(C), A, B) ≈ check1
+    check = mul!(Array(C), Array(A), Array(B))
+    @test mul!(copy(C), A, B) ≈ check
+
+    check = mul!(Array(C), Array(A)', Array(B))
+    @test mul!(copy(C), A', B) ≈ check
+
+    check = mul!(Array(C), Array(A), Array(B)')
+    @test mul!(copy(C), A, B') ≈ check
+
+    check = mul!(Array(C), Array(A)', Array(B)')
+    @test mul!(copy(C), A', B') ≈ check
 
     α = rand(rng, T)
     β = rand(rng, T)
-    check2 = mul!(Array(C), Array(A), Array(B), α, β)
-    @test mul!(copy(C), A, B, α, β) ≈ check2
+    check = mul!(Array(C), Array(A), Array(B), α, β)
+    @test mul!(copy(C), A, B, α, β) ≈ check
   end
 
   # test empty matrix