Fix issue with complex adjoint views (#14)

* Fix issue with complex adjoint views

* Introduce ConjPtr

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ArrayLayouts"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.2.3"
+version = "0.2.4"
 
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"

src/ArrayLayouts.jl

@@ -27,7 +27,7 @@ import Base: AbstractArray, AbstractMatrix, AbstractVector,
       similar, @_gc_preserve_end, @_gc_preserve_begin,
       @nexprs, @ncall, @ntuple, tuple_type_tail,
       all, any, isbitsunion, issubset, replace_in_print_matrix, replace_with_centered_mark,
-      strides, unsafe_convert
+      strides, unsafe_convert, first_index
 
 import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcasted,
                         combine_eltypes, DefaultArrayStyle, instantiate, materialize,
@@ -71,8 +71,24 @@ _transpose_strides(a,b) = (b,a)
 strides(A::Adjoint) = _transpose_strides(strides(parent(A))...)
 strides(A::Transpose) = _transpose_strides(strides(parent(A))...)
 
+"""
+    ConjPtr{T}
+
+represents that the entry is the complex-conjugate of the pointed to entry.
+"""
+struct ConjPtr{T} 
+    ptr::Ptr{T}
+end
+
 unsafe_convert(::Type{Ptr{T}}, A::Adjoint{<:Real}) where T<:Real = unsafe_convert(Ptr{T}, parent(A))
 unsafe_convert(::Type{Ptr{T}}, A::Transpose) where T = unsafe_convert(Ptr{T}, parent(A))
+# work-around issue with complex conjugation of pointer
+unsafe_convert(::Type{Ptr{T}}, Ac::Adjoint{<:Complex}) where T<:Complex = unsafe_convert(ConjPtr{T}, parent(Ac))
+unsafe_convert(::Type{ConjPtr{T}}, Ac::Adjoint{<:Complex}) where T<:Complex = unsafe_convert(Ptr{T}, parent(Ac))    
+function unsafe_convert(::Type{ConjPtr{T}}, V::SubArray{T,2}) where {T,N,P}
+    kr, jr = parentindices(V)
+    unsafe_convert(Ptr{T}, view(parent(V)', jr, kr))
+end
 
 include("memorylayout.jl")
 include("muladd.jl")

src/muladd.jl

@@ -257,64 +257,65 @@ end
     end
 end
 
+# work around pointer issues
 @inline materialize!(M::BlasMatMulVecAdd{<:AbstractColumnMajor,<:AbstractStridedLayout,<:AbstractStridedLayout}) =
     _gemv!('N', M.α, M.A, M.B, M.β, M.C)
 @inline materialize!(M::BlasMatMulVecAdd{<:AbstractRowMajor,<:AbstractStridedLayout,<:AbstractStridedLayout}) =
     _gemv!('T', M.α, transpose(M.A), M.B, M.β, M.C)
 @inline materialize!(M::BlasMatMulVecAdd{<:ConjLayout{<:AbstractRowMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasComplex}) =
-    _gemv!('C', M.α, M.A', M.B, M.β, M.C)
+    _gemv!('C', M.α, adjoint(M.A), M.B, M.β, M.C)
 
 @inline materialize!(M::BlasVecMulMatAdd{<:AbstractColumnMajor,<:AbstractColumnMajor,<:AbstractColumnMajor}) =
     _gemm!('N', 'N', M.α, M.A, M.B, M.β, M.C)
 @inline materialize!(M::BlasVecMulMatAdd{<:AbstractColumnMajor,<:AbstractRowMajor,<:AbstractColumnMajor}) =
     _gemm!('N', 'T', M.α, M.A, transpose(M.B), M.β, M.C)
 @inline materialize!(M::BlasVecMulMatAdd{<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
-    _gemm!('N', 'C', M.α, M.A, M.B', M.β, M.C)
+    _gemm!('N', 'C', M.α, M.A, adjoint(M.B), M.β, M.C)
 
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractColumnMajor,<:AbstractColumnMajor}) =
     _gemm!('N', 'N', M.α, M.A, M.B, M.β, M.C)
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractRowMajor,<:AbstractColumnMajor}) =
     _gemm!('N', 'T', M.α, M.A, transpose(M.B), M.β, M.C)
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
-    _gemm!('N', 'C', M.α, M.A, M.B', M.β, M.C)
+    _gemm!('N', 'C', M.α, M.A, adjoint(M.B), M.β, M.C)
 
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:AbstractColumnMajor,<:AbstractColumnMajor}) =
     _gemm!('T', 'N', M.α, transpose(M.A), M.B, M.β, M.C)
 @inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:AbstractColumnMajor,<:BlasComplex}) =
-    _gemm!('C', 'N', M.α, M.A', M.B, M.β, M.C)
+    _gemm!('C', 'N', M.α, adjoint(M.A), M.B, M.β, M.C)
 
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:AbstractRowMajor,<:AbstractColumnMajor}) =
     _gemm!('T', 'T', M.α, transpose(M.A), transpose(M.B), M.β, M.C)
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
-    _gemm!('T', 'C', M.α, transpose(M.A), M.B', M.β, M.C)
+    _gemm!('T', 'C', M.α, transpose(M.A), adjoint(M.B), M.β, M.C)
 
 @inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:AbstractRowMajor,<:AbstractColumnMajor,<:BlasComplex}) =
-    _gemm!('C', 'T', M.α, M.A', M.B', M.β, M.C)
+    _gemm!('C', 'T', M.α, adjoint(M.A), transpose(M.B), M.β, M.C)
 @inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
-    _gemm!('C', 'C', M.α, M.A', M.B', M.β, M.C)
+    _gemm!('C', 'C', M.α, adjoint(M.A), adjoint(M.B), M.β, M.C)
 
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractColumnMajor,<:AbstractRowMajor}) =
     _gemm!('T', 'T', M.α, M.B, M.A, M.β, transpose(M.C))
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
-    _gemm!('C', 'C', M.α, M.B, M.A, M.β, M.C')
+    _gemm!('C', 'C', M.α, M.B, M.A, M.β, adjoint(M.C))
 
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractRowMajor,<:AbstractRowMajor}) =
     _gemm!('N', 'T', M.α, transpose(M.B), M.A, M.β, transpose(M.C))
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractRowMajor,<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
     _gemm!('N', 'T', M.α, transpose(M.B), M.A, M.β, M.C')
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
-    _gemm!('N', 'C', M.α, M.B', M.A, M.β, M.C')
+    _gemm!('N', 'C', M.α, adjoint(M.B), M.A, M.β, adjoint(M.C))
 
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:AbstractColumnMajor,<:AbstractRowMajor}) =
     _gemm!('T', 'N', M.α, M.B, transpose(M.A), M.β, transpose(M.C))
 @inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
-    _gemm!('C', 'N', M.α, M.B, M.A', M.β, M.C')
+    _gemm!('C', 'N', M.α, M.B, adjoint(M.A), M.β, adjoint(M.C))
 
 
 @inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:AbstractRowMajor,<:AbstractRowMajor}) =
     _gemm!('N', 'N', M.α, transpose(M.B), transpose(M.A), M.β, transpose(M.C))
 @inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
-    _gemm!('N', 'N', M.α, M.B', M.A', M.β, M.C')
+    _gemm!('N', 'N', M.α, adjoint(M.B), adjoint(M.A), M.β, adjoint(M.C))
 
 
 ###

test/test_muladd.jl

@@ -168,6 +168,12 @@ end
             c .= MulAdd(3one(T),Ac,b,2one(T),c)
             @test all(c .=== BLAS.gemv!(trans, 3one(T), A, b, 2one(T), copy(b)))
         end
+
+        C = randn(6,6) + im*randn(6,6)
+        V = view(C', 2:3, 3:4)
+        c = randn(2) + im*randn(2)
+        @test all(muladd!(1.0+0im,V,c,0.0+0im,similar(c,2)) .=== BLAS.gemv!('C', 1.0+0im, Matrix(V'), c, 0.0+0im, similar(c,2)))
+        @test all(muladd!(1.0+0im,V',c,0.0+0im,similar(c,2)) .=== BLAS.gemv!('N', 1.0+0im, Matrix(V'), c, 0.0+0im, similar(c,2)))
     end
 
     @testset "gemm adjtrans" begin