Make transposes of StridedArrays strided

[mbauman]

but only in cases where the transpose is actually stored in memory.

[dlfivefifty]

Would it be at all beneficial to use reverse(strides(parent(A))) ?

[mbauman]

Sure enough:

julia> t = transpose(rand(3,4));

julia> s1(A) = (stride(A.parent, 2), stride(A.parent, 1))
       s2(A) = reverse(strides(A.parent))
s2 (generic function with 1 method)

julia> @btime s1($t)
  7.056 ns (0 allocations: 0 bytes)
(3, 1)

julia> @btime s2($t)
  2.032 ns (0 allocations: 0 bytes)
(3, 1)

[mbauman]

Oh wait, that won't quite do it because of transposed vectors… should have actually tested it more thoroughly locally.

[dlfivefifty]

_adjoint_strides((a,)) = (a,1)
_adjoint_strides((a,b)) = (b,a)
strides(A::Adoint) = _adjoint_strides(strides(parent(A)))

[mbauman]

The answer for strided vectors is neither reverse(strides(A)) nor (stride(A,1), 1) nor (length(A), stride(A, 1)). Given that all three of us got this optimization wrong at first glance, I'm heavily inclined to just ask the vector itself for its stride(A, 2) instead of trying to re-implement that computation ourselves as a micro-optimization.

(I believe the correct answer is s = stride(A, 1); return (length(A)*s, s), but given that you should only ever multiply this number by 0 the point is somewhat moot… although you might get yourself in trouble with some error-checking routines that assert the strides are non-overlapping).

[mbauman]

Let's re-run CI and get this in.

[mcabbott]

Perhaps related, pointer(transpose(ones(2,2)) gives an error right now (on Julia 1.4), should this also work?

[mbauman]

Yes, I think we could add pointer(::Transpose) (and it make sense to do so in this PR) but I don't think we should add pointer(::Adjoint) as the latter elements aren't actually stored on disk.

[mcabbott]

Shouldn't it follow strides and work for an Adjoint of real numbers?

[mbauman]

Ah, yes, I had forgotten I did that. That's perfect.

[dlfivefifty]

This should not restrict to <: StridedVector. The strided array interface means other array types should be supported.

[mbauman]

Yes, I fully agree, but this is intended as a stopgap and not a full solution.

[dlfivefifty]

FYI In ArrayLayouts.jl I came up with the following way to deal with the complex adjoint pointer issue: introduce ConjPtr to represent that the pointer is to the conjugate of the element.

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

[mcabbott]

Bump, for 1.5?

[mbauman]

So was this discussed on the triage call? I really don't think it should be all that controversial — the only controversial thing is how minimally I've scoped it in an effort to keep in uncontroversial.

[StefanKarpinski]

Your call, @mbauman.

[mbauman]

Let's do this then. Part of the reason it's scoped so narrowly is to ensure we have space to move towards traits or whatever the bigger solution turns out to be.