Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make transposes of StridedArrays strided #29135

Merged
merged 5 commits into from
May 4, 2020
Merged

Make transposes of StridedArrays strided #29135

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
1032431
Make transposes of StridedArrays strided
mbauman Sep 11, 2018
57cd458
Use reverse instead
mbauman Sep 14, 2018
b6ff2c9
Split definitions for vectors and matrices
mbauman Sep 14, 2018
a9a8af7
Similarly add pointer support
mbauman Apr 1, 2020
9819783
Use unsafe_convert instead of pointer
mbauman Apr 28, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions stdlib/LinearAlgebra/src/adjtrans.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -186,6 +186,16 @@ IndexStyle(::Type{<:AdjOrTransAbsMat}) = IndexCartesian()
convert(::Type{Adjoint{T,S}}, A::Adjoint) where {T,S} = Adjoint{T,S}(convert(S, A.parent))
convert(::Type{Transpose{T,S}}, A::Transpose) where {T,S} = Transpose{T,S}(convert(S, A.parent))

# Strides and pointer for transposed strided arrays — but only if the elements are actually stored in memory
Base.strides(A::Adjoint{<:Real, <:StridedVector}) = (stride(A.parent, 2), stride(A.parent, 1))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Base.strides(A::Transpose{<:Any, <:StridedVector}) = (stride(A.parent, 2), stride(A.parent, 1))
# For matrices it's slightly faster to use reverse and avoid calling stride twice
Base.strides(A::Adjoint{<:Real, <:StridedMatrix}) = reverse(strides(A.parent))
Base.strides(A::Transpose{<:Any, <:StridedMatrix}) = reverse(strides(A.parent))

Base.unsafe_convert(::Type{Ptr{T}}, A::Adjoint{<:Real, <:StridedVecOrMat}) where {T} = Base.unsafe_convert(Ptr{T}, A.parent)
Base.unsafe_convert(::Type{Ptr{T}}, A::Transpose{<:Any, <:StridedVecOrMat}) where {T} = Base.unsafe_convert(Ptr{T}, A.parent)

# for vectors, the semantics of the wrapped and unwrapped types differ
# so attempt to maintain both the parent and wrapper type insofar as possible
similar(A::AdjOrTransAbsVec) = wrapperop(A)(similar(A.parent))
Expand Down
25 changes: 25 additions & 0 deletions stdlib/LinearAlgebra/test/adjtrans.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -483,6 +483,31 @@ end
"$t of "*sprint((io, t) -> show(io, MIME"text/plain"(), t), parent(Fop))
end

@testset "strided transposes" begin
for t in (Adjoint, Transpose)
@test strides(t(rand(3))) == (3, 1)
@test strides(t(rand(3,2))) == (3, 1)
@test strides(t(view(rand(3, 2), :))) == (6, 1)
@test strides(t(view(rand(3, 2), :, 1:2))) == (3, 1)

A = rand(3)
@test pointer(t(A)) === pointer(A)
B = rand(3,1)
@test pointer(t(B)) === pointer(B)
end
@test_throws MethodError strides(Adjoint(rand(3) .+ rand(3).*im))
@test_throws MethodError strides(Adjoint(rand(3, 2) .+ rand(3, 2).*im))
@test strides(Transpose(rand(3) .+ rand(3).*im)) == (3, 1)
@test strides(Transpose(rand(3, 2) .+ rand(3, 2).*im)) == (3, 1)

C = rand(3) .+ rand(3).*im
@test_throws ErrorException pointer(Adjoint(C))
@test pointer(Transpose(C)) === pointer(C)
D = rand(3,2) .+ rand(3,2).*im
@test_throws ErrorException pointer(Adjoint(D))
@test pointer(Transpose(D)) === pointer(D)
end

const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
isdefined(Main, :OffsetArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "OffsetArrays.jl"))
using .Main.OffsetArrays
Expand Down