From 18b378d40e2d3739053ad4c936ca2fab697e7a0a Mon Sep 17 00:00:00 2001
From: Valentin Churavy <v.churavy@gmail.com>
Date: Tue, 14 Apr 2020 08:32:51 -0700
Subject: [PATCH 1/9] Accept primitive types and multiple array types

Changes:
- Accepts primitive types (e.g. types with no fields).
  This allows broadcast to produce SoA arrays starting
  from something like `SoA{Float32, 2}`.
- Supports multiple storage array types. Most notably
  `CuArrays` and `CuDeviceArray`.
---
 LICENSE.md             |   1 +
 src/StructsOfArrays.jl | 184 +++++++++++++++++++++++++++++++----------
 2 files changed, 141 insertions(+), 44 deletions(-)

diff --git a/LICENSE.md b/LICENSE.md
index a3451f8..c2fc6d1 100644
--- a/LICENSE.md
+++ b/LICENSE.md
@@ -1,6 +1,7 @@
 The StructsOfArrays.jl package is licensed under the MIT "Expat" License:
 
 > Copyright (c) 2015: Simon Kornblith.
+> Copyright (c) 2018-2019: Valentin Churavy, and other contributors
 >
 > Permission is hereby granted, free of charge, to any person obtaining
 > a copy of this software and associated documentation files (the
diff --git a/src/StructsOfArrays.jl b/src/StructsOfArrays.jl
index eb4226a..048cd46 100644
--- a/src/StructsOfArrays.jl
+++ b/src/StructsOfArrays.jl
@@ -1,16 +1,49 @@
 module StructsOfArrays
 
-using Compat
-
 export StructOfArrays
 
-immutable StructOfArrays{T,N,U<:Tuple} <: AbstractArray{T,N}
+using Adapt
+
+"""
+    StructOfArrays{T, N, AT, U} <: AbstractArray{T,N}
+
+Transparent Struct of Arrays transformation.
+Given an array like `Array{ComplexF64, 2}`, we want to represent it as:
+
+```
+struct ComplexArray{N}
+    re::Array{Float64}
+    im::Array{Float64}
+end
+
+getindex(A::ComplexArray, i) = ComplexF64(A.re[i], A.im[i])
+```
+
+Since the memory-access patterns of `ComplexArray` are much more SIMD friendly than
+`Array{ComplexF64}`. This version of `StructOfArrays` also has a notion of underlying
+storage and works both on CPU arrays and CuArrays. In most cases single field access of
+the struct should be fast as well since Julia and LLVM both do DCE.
+"""
+struct StructOfArrays{T,N, AT<:AbstractArray{T, N},U<:Tuple} <: AbstractArray{T,N}
     arrays::U
 end
 
+# Storage types of StructOfArrays need to implement this
+_type_with_eltype(::Type{<:Array}, T, N) = Array{T, N}
+_type(::Type{<:Array}) = Array
+
+function Adapt.adapt_structure(to, x::StructOfArrays{T, N}) where {T, N}
+    arrays = map(A -> adapt(to, A), x.arrays)
+    TT = typeof(arrays)
+    AT = _type_with_eltype(to, T, N)
+    StructOfArrays{T, N, AT, TT}(arrays)
+end
+
 function gather_eltypes(T, visited = Set{Type}())
-    (!isleaftype(T) || T.mutable) && throw(ArgumentError("can only create an StructOfArrays of leaf type immutables"))
-    isempty(T.types) && throw(ArgumentError("cannot create an StructOfArrays of an empty or bitstype"))
+    (!isconcretetype(T) || T.mutable) && throw(ArgumentError("can only create an StructOfArrays of leaf type immutables"))
+    if isempty(T.types)
+        return Type[T]
+    end
     types = Type[]
     push!(visited, T)
     for S in T.types
@@ -25,59 +58,68 @@ function gather_eltypes(T, visited = Set{Type}())
     types
 end
 
-@generated function StructOfArrays{T}(::Type{T}, dims::Integer...)
-    N = length(dims)
-    types = gather_eltypes(T)
-    arrtuple = Tuple{[Array{S,N} for S in types]...}
-    :(StructOfArrays{T,$N,$arrtuple}(($([:(Array{$(S)}(dims)) for S in types]...),)))
-end
-StructOfArrays(T::Type, dims::Tuple{Vararg{Integer}}) = StructOfArrays(T, dims...)
+@generated function StructOfArrays(::Type{T}, ::Type{ArrayT}, dims::Integer...) where {T, ArrayT<:AbstractArray}
+    N         = length(dims)
+    pArrayT   = _type_with_eltype(ArrayT, T, N)
+    types     = gather_eltypes(T)
+    arrtypes  = map(t->_type_with_eltype(ArrayT, t, N), types)
+    arrtuple  = Tuple{arrtypes...}
 
-@compat Base.IndexStyle{T<:StructOfArrays}(::Type{T}) = IndexLinear()
-
-@generated function Base.similar{T}(A::StructOfArrays, ::Type{T}, dims::Dims)
-    if isbits(T) && length(T.types) > 1
-        :(StructOfArrays(T, dims))
-    else
-        :(Array{T}(dims))
-    end
+    :(StructOfArrays{T,$N,$(pArrayT),$arrtuple}(($([:($(arrtypes[i])(undef,dims)) for i = 1:length(types)]...),)))
 end
+StructOfArrays(T::Type, AT::Type, dims::Tuple{Vararg{Integer}}) = StructOfArrays(T, AT, dims...)
 
-Base.convert{T,S,N}(::Type{StructOfArrays{T,N}}, A::AbstractArray{S,N}) =
-    copy!(StructOfArrays(T, size(A)), A)
-Base.convert{T,S,N}(::Type{StructOfArrays{T}}, A::AbstractArray{S,N}) =
-    convert(StructOfArrays{T,N}, A)
-Base.convert{T,N}(::Type{StructOfArrays}, A::AbstractArray{T,N}) =
-    convert(StructOfArrays{T,N}, A)
+Base.size(A::StructOfArrays) = size(@inbounds(A.arrays[1]))
+Base.size(A::StructOfArrays, i) = size(@inbounds(A.arrays[1]), i)
+
+Base.show(io::IO, a::StructOfArrays{T,N,A}) where {T,N,A} = print(io, "$(length(a))-element SoA{$T,$N,$A}")
 
-Base.size(A::StructOfArrays) = size(A.arrays[1])
-Base.size(A::StructOfArrays, d) = size(A.arrays[1], d)
+Base.print_array(::IO, ::StructOfArrays) = nothing
 
-function generate_getindex(T, arraynum)
+function generate_getindex(T, getindex, arraynum)
     members = Expr[]
     for S in T.types
-        sizeof(S) == 0 && push!(members, :($(S())))
-        if isempty(S.types)
-            push!(members, :(A.arrays[$arraynum][i...]))
+        if sizeof(S) == 0
+            if S <: Tuple
+                exprs2 = Expr[]
+                for S2 in S.parameters
+                    push!(exprs2, :($(S2)()))
+                end
+                push!(members, Expr(:tuple, exprs2))
+            else
+                push!(members, :($(S)()))
+            end
+        elseif isempty(S.types)
+            push!(members, :($(getindex)(A.arrays[$arraynum], i...)))
             arraynum += 1
         else
-            member, arraynum = generate_getindex(S, arraynum)
+            member, arraynum = generate_getindex(S, getindex, arraynum)
             push!(members, member)
         end
     end
     Expr(:new, T, members...), arraynum
 end
 
-@generated function Base.getindex{T}(A::StructOfArrays{T}, i::Integer...)
-    strct, _ = generate_getindex(T, 1)
-    Expr(:block, Expr(:meta, :inline), strct)
+@generated function Base.getindex(A::StructOfArrays{T}, i::Integer...) where T
+    exprs = Any[Expr(:meta, :inline), Expr(:meta, :propagate_inbounds)]
+    if isempty(T.types)
+        push!(exprs, :(return A.arrays[1][i...]))
+    else
+        strct, _ = generate_getindex(T, Base.getindex, 1)
+        push!(exprs, strct)
+    end
+    quote
+        $(exprs...)
+    end
 end
 
 function generate_setindex(T, x, arraynum)
     s = gensym()
     exprs = Expr[:($s = $x)]
     for (el,S) in enumerate(T.types)
-        sizeof(S) == 0 && push!(members, :($(S())))
+        if sizeof(S) == 0
+            continue
+        end
         if isempty(S.types)
             push!(exprs, :(A.arrays[$arraynum][i...] = getfield($s, $el)))
             arraynum += 1
@@ -89,13 +131,67 @@ function generate_setindex(T, x, arraynum)
     exprs, arraynum
 end
 
-@generated function Base.setindex!{T}(A::StructOfArrays{T}, x, i::Integer...)
-    exprs = Expr(:block, generate_setindex(T, :x, 1)[1]...)
+@generated function Base.setindex!(A::StructOfArrays{T}, x, i::Integer...) where {T}
+    exprs = Any[Expr(:meta, :inline), Expr(:meta, :propagate_inbounds)]
+    push!(exprs, :(v = convert(T, x)))
+    if isempty(T.types)
+        push!(exprs, :(A.arrays[1][i...] = v))
+    else
+        append!(exprs, generate_setindex(T, :v, 1)[1])
+    end
+    push!(exprs, :(return x))
     quote
-        $(Expr(:meta, :inline))
-        v = convert(T, x)
-        $exprs
-        x
+        $(exprs...)
+    end
+end
+
+# Base.IndexStyle(::Type{<:StructOfArrays{T,N,A}}) where {T,N,A<:AbstractArray} =  Base.IndexStyle(A)
+# Base.BroadcastStyle(::Type{<:StructOfArrays{T,N,A}}) where {T,N,A<:AbstractArray} = Broadcast.ArrayStyle{StructOfArrays{T,N,A}}()
+
+# function Base.similar(A::StructOfArrays{T1,N,AT}, ::Type{T}, dims::Dims) where {T1,N,AT,T}
+#     StructOfArrays(T, AT, dims)
+# end
+
+# function Base.broadcast_similar(f, ::Broadcast.ArrayStyle{StructOfArrays{T1,N,A}}, ::Type{T}, inds, As...) where {T1,N,A,T}
+#     StructOfArrays(T, A, Base.to_shape(inds))
+# end
+
+function Base.convert(::Type{<:StructOfArrays{T,N,AT}}, A::StructOfArrays{T,N}) where {T,N,AT<:AbstractArray{T,N}}
+    if AT <: StructOfArrays
+        error("Can't embed a SoA array in a SoA array")
     end
+    arrays = map(a->convert(_type(AT), a), A.arrays)
+    tt = typeof(arrays)
+    StructOfArrays{T, N, AT, tt}(arrays)
+end
+
+function Base.convert(::Type{<:StructOfArrays{T,N,AT}}, A::StructOfArrays{S,N,BT}) where {T,N,AT,S,BT}
+    BT != AT && AT<:CuArray && error("Can't convert from $BT to $AT with different eltypes")
+    copyto!(StructOfArrays(T, _type_with_eltype(AT, T, N), size(A)), A)
+end
+
+function Base.convert(::Type{<:StructOfArrays{T,N}}, A::AbstractArray) where {T,N}
+    @assert !(A isa StructOfArrays)
+    copyto!(StructOfArrays(T, _type_with_eltype(typeof(A), T, N), size(A)), A)
+end
+
+Base.convert(::Type{<:StructOfArrays{T}}, A::AbstractArray{S,N}) where {T,S,N} =
+    convert(StructOfArrays{T,N}, A)
+
+Base.convert(::Type{<:StructOfArrays}, A::AbstractArray{T,N}) where {T,N} =
+    convert(StructOfArrays{T,N}, A)
+
+function Base.one(x::StructOfArrays{T}) where {T}
+    # FIXME: probably not the best way, ie. don't use Base._one
+    convert(StructOfArrays, Base._one(one(T), x))
+end
+
+function Base.Array{T,N}(src::StructOfArrays{T,N}) where {T,N}
+    A = convert(StructOfArrays{T,N,Array{T,N}}, src)
+    dst = Array{T,N}(uninitialized, size(src))
+    copyto!(dst, A)
+    return dst
+end
+Array(src::StructOfArrays{T,N}) where {T,N} = Array{T,N}(src)
+
 end
-end # module

From db6366141239fafe182b78b5152143a4b46e6c64 Mon Sep 17 00:00:00 2001
From: Lucas C Wilcox <lucas@swirlee.com>
Date: Tue, 14 Apr 2020 10:04:17 -0700
Subject: [PATCH 2/9] Make this a julia package

In this process the broadcasting was updated to work and most of the
tests were updated to run.
---
 .gitignore             |  5 ++-
 .travis.yml            | 16 +++----
 LICENSE                | 20 +++++++++
 LICENSE.md             | 23 ----------
 Project.toml           | 11 +++++
 REQUIRE                |  2 -
 appveyor.yml           | 37 ----------------
 src/StructsOfArrays.jl | 16 +++----
 test/Project.toml      |  4 ++
 test/REQUIRE           |  1 -
 test/runtests.jl       | 95 +++++++++++++++++++++++-------------------
 11 files changed, 108 insertions(+), 122 deletions(-)
 create mode 100644 LICENSE
 delete mode 100644 LICENSE.md
 create mode 100644 Project.toml
 delete mode 100644 REQUIRE
 delete mode 100644 appveyor.yml
 create mode 100644 test/Project.toml
 delete mode 100644 test/REQUIRE

diff --git a/.gitignore b/.gitignore
index 42f07ea..2753f7e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,3 @@
-*.jl.cov
-*.jl.mem
+.DS_Store
+/Manifest.toml
+/test/Manifest.toml
diff --git a/.travis.yml b/.travis.yml
index 1ce59d4..29e6ea3 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,15 +1,17 @@
+## Documentation: http://docs.travis-ci.com/user/languages/julia/
 language: julia
 os:
   - linux
   - osx
+  - windows
 julia:
-  - 0.4
-  - 0.5
+  - 1.3
+  - 1.4
   - nightly
 notifications:
   email: false
-script:
- - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
- - julia --check-bounds=yes -e 'Pkg.clone(pwd()); Pkg.build("StructsOfArrays"); Pkg.test("StructsOfArrays"; coverage=true)'
-after_success:
-  - julia -e 'cd(Pkg.dir("StructsOfArrays")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())'
+git:
+  depth: 99999999
+jobs:
+  allow_failures:
+    - julia: nightly
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..29eadaa
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,20 @@
+Copyright (c) 2015: Simon Kornblith.
+Copyright (c) 2018-2019: Valentin Churavy, and other contributors
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/LICENSE.md b/LICENSE.md
deleted file mode 100644
index c2fc6d1..0000000
--- a/LICENSE.md
+++ /dev/null
@@ -1,23 +0,0 @@
-The StructsOfArrays.jl package is licensed under the MIT "Expat" License:
-
-> Copyright (c) 2015: Simon Kornblith.
-> Copyright (c) 2018-2019: Valentin Churavy, and other contributors
->
-> Permission is hereby granted, free of charge, to any person obtaining
-> a copy of this software and associated documentation files (the
-> "Software"), to deal in the Software without restriction, including
-> without limitation the rights to use, copy, modify, merge, publish,
-> distribute, sublicense, and/or sell copies of the Software, and to
-> permit persons to whom the Software is furnished to do so, subject to
-> the following conditions:
->
-> The above copyright notice and this permission notice shall be
-> included in all copies or substantial portions of the Software.
->
-> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-> EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-> MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-> IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
-> CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
-> TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
-> SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
diff --git a/Project.toml b/Project.toml
new file mode 100644
index 0000000..032ff7c
--- /dev/null
+++ b/Project.toml
@@ -0,0 +1,11 @@
+name = "StructsOfArrays"
+uuid = "ad108753-388a-4a68-a626-efa2fa9e4278"
+authors = ["Simon Kornblith <simon@simonster.com>", "Valentin Churavy <v.churavy@gmail.com>", "Other Contributors"]
+version = "0.0.3"
+
+[deps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+
+[compat]
+Adapt = "~1"
+julia = "≥ 1.3"
diff --git a/REQUIRE b/REQUIRE
deleted file mode 100644
index 33ed633..0000000
--- a/REQUIRE
+++ /dev/null
@@ -1,2 +0,0 @@
-julia 0.4
-Compat 0.19
diff --git a/appveyor.yml b/appveyor.yml
deleted file mode 100644
index 979983a..0000000
--- a/appveyor.yml
+++ /dev/null
@@ -1,37 +0,0 @@
-environment:
-  matrix:
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.4/julia-0.4-latest-win32.exe"
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.4/julia-0.4-latest-win64.exe"
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.5/julia-0.5-latest-win32.exe"
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.5/julia-0.5-latest-win64.exe"
-  - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
-  - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"
-
-branches:
-  only:
-    - master
-    - /release-.*/
-
-notifications:
-  - provider: Email
-    on_build_success: false
-    on_build_failure: false
-    on_build_status_changed: false
-
-install:
-  - ps: "[System.Net.ServicePointManager]::SecurityProtocol = [System.Net.SecurityProtocolType]::Tls12"
-# Download most recent Julia Windows binary
-  - ps: (new-object net.webclient).DownloadFile(
-        $env:JULIA_URL,
-        "C:\projects\julia-binary.exe")
-# Run installer silently, output to C:\projects\julia
-  - C:\projects\julia-binary.exe /S /D=C:\projects\julia
-
-build_script:
-# Need to convert from shallow to complete for Pkg.clone to work
-  - IF EXIST .git\shallow (git fetch --unshallow)
-  - C:\projects\julia\bin\julia -e "versioninfo();
-      Pkg.clone(pwd(), \"StructsOfArrays\"); Pkg.build(\"StructsOfArrays\")"
-
-test_script:
-  - C:\projects\julia\bin\julia --check-bounds=yes -e "Pkg.test(\"StructsOfArrays\")"
diff --git a/src/StructsOfArrays.jl b/src/StructsOfArrays.jl
index 048cd46..6ccfaf8 100644
--- a/src/StructsOfArrays.jl
+++ b/src/StructsOfArrays.jl
@@ -145,16 +145,16 @@ end
     end
 end
 
-# Base.IndexStyle(::Type{<:StructOfArrays{T,N,A}}) where {T,N,A<:AbstractArray} =  Base.IndexStyle(A)
-# Base.BroadcastStyle(::Type{<:StructOfArrays{T,N,A}}) where {T,N,A<:AbstractArray} = Broadcast.ArrayStyle{StructOfArrays{T,N,A}}()
+Base.IndexStyle(::Type{<:StructOfArrays{T,N,A}}) where {T,N,A<:AbstractArray} =  Base.IndexStyle(A)
+Base.BroadcastStyle(::Type{<:StructOfArrays{T,N,A}}) where {T,N,A<:AbstractArray} = Broadcast.ArrayStyle{StructOfArrays{T,N,A}}()
 
-# function Base.similar(A::StructOfArrays{T1,N,AT}, ::Type{T}, dims::Dims) where {T1,N,AT,T}
-#     StructOfArrays(T, AT, dims)
-# end
+function Base.similar(A::StructOfArrays{T1,N,AT}, ::Type{T}, dims::Dims) where {T1,N,AT,T}
+    StructOfArrays(T, AT, dims)
+end
 
-# function Base.broadcast_similar(f, ::Broadcast.ArrayStyle{StructOfArrays{T1,N,A}}, ::Type{T}, inds, As...) where {T1,N,A,T}
-#     StructOfArrays(T, A, Base.to_shape(inds))
-# end
+function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{StructOfArrays{T1,N,AT}}}, ::Type{T}) where {T1,N,AT,T}
+    StructOfArrays(T, AT, Base.to_shape(axes(bc)))
+end
 
 function Base.convert(::Type{<:StructOfArrays{T,N,AT}}, A::StructOfArrays{T,N}) where {T,N,AT<:AbstractArray{T,N}}
     if AT <: StructOfArrays
diff --git a/test/Project.toml b/test/Project.toml
new file mode 100644
index 0000000..5342819
--- /dev/null
+++ b/test/Project.toml
@@ -0,0 +1,4 @@
+[deps]
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StructsOfArrays = "ad108753-388a-4a68-a626-efa2fa9e4278"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/REQUIRE b/test/REQUIRE
deleted file mode 100644
index fbb6170..0000000
--- a/test/REQUIRE
+++ /dev/null
@@ -1 +0,0 @@
-Compat 0.9.1
diff --git a/test/runtests.jl b/test/runtests.jl
index 98b30bc..91aa5d6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,44 +1,55 @@
 using StructsOfArrays
-using Base.Test
-using Compat
-
-regular = @compat complex.(randn(10000), randn(10000))
-soa = convert(StructOfArrays, regular)
-@test regular == soa
-@test sum(regular) ≈ sum(soa)
-
-soa64 = convert(StructOfArrays{Complex64}, regular)
-@test convert(Array{Complex64}, regular) == soa64
-
-sim = similar(soa)
-@test typeof(sim) == typeof(soa)
-@test size(sim) == size(soa)
-
-regular = @compat complex.(randn(10, 5), randn(10, 5))
-soa = convert(StructOfArrays, regular)
-for i = 1:10, j = 1:5
-    @test regular[i, j] == soa[i, j]
+using Test
+using Random
+
+@testset "StructsOfArrays.jl" begin
+    @testset "constructor" begin
+        regular = rand(MersenneTwister(0), ComplexF64, 10000)
+        soa = convert(StructOfArrays, regular)
+        @test regular == soa
+        @test sum(regular) ≈ sum(soa)
+
+        soa64 = convert(StructOfArrays{ComplexF64}, regular)
+        @test convert(Array{ComplexF64}, regular) == soa64
+
+        sim = similar(soa)
+        @test typeof(sim) == typeof(soa)
+        @test size(sim) == size(soa)
+
+        regular = randn(MersenneTwister(0), ComplexF64, 10, 5)
+        soa = convert(StructOfArrays, regular)
+        for i = 1:10, j = 1:5
+            @test regular[i, j] == soa[i, j]
+        end
+        @test size(soa, 1) == 10
+        @test size(soa, 2) == 5
+    end
+
+    @testset "similar" begin
+        struct OneField
+            x::Int
+        end
+
+        small = StructOfArrays(ComplexF64, Array, 2)
+        @test typeof(small) <: AbstractArray{Complex{T}} where T
+        @test typeof(similar(small, ComplexF64)) <: AbstractArray{Complex{Float64}}
+        @test typeof(similar(small, Int)) <: AbstractArray{Int}
+        @test typeof(similar(small, OneField)) <: AbstractArray{OneField}
+        @test typeof(similar(small, ComplexF64)) <: StructOfArrays
+    end
+
+    @testset "recursive structs" begin
+        struct OneField
+            x::Int
+        end
+
+        struct TwoField
+            one::OneField
+            two::OneField
+        end
+
+        small = StructOfArrays(TwoField, Array, 2, 2)
+        small[1,1] = TwoField(OneField(1), OneField(2))
+        @test small[1,1] == TwoField(OneField(1), OneField(2))
+    end
 end
-@test size(soa, 1) == 10
-@test size(soa, 2) == 5
-
-immutable OneField
-    x::Int
-end
-
-small = StructOfArrays(Complex64, 2)
-@test typeof(similar(small, Complex)) === Vector{Complex}
-@test typeof(similar(small, Int)) === Vector{Int}
-@test typeof(similar(small, SubString)) === Vector{SubString}
-@test typeof(similar(small, OneField)) === Vector{OneField}
-@test typeof(similar(small, Complex128)) <: StructOfArrays
-
-# Recursive structs
-immutable TwoField
-    one::OneField
-    two::OneField
-end
-
-small = StructOfArrays(TwoField, 2, 2)
-small[1,1] = TwoField(OneField(1), OneField(2))
-@test small[1,1] == TwoField(OneField(1), OneField(2))

From 8f2e83a94755793c886a03781f310a2f9ab21943 Mon Sep 17 00:00:00 2001
From: Lucas C Wilcox <lucas@swirlee.com>
Date: Tue, 14 Apr 2020 13:07:40 -0700
Subject: [PATCH 3/9] Add basic test with eltype of StaticArray

---
 test/Project.toml |  1 +
 test/runtests.jl  | 11 +++++++++++
 2 files changed, 12 insertions(+)

diff --git a/test/Project.toml b/test/Project.toml
index 5342819..e91b19a 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,4 +1,5 @@
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StructsOfArrays = "ad108753-388a-4a68-a626-efa2fa9e4278"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/runtests.jl b/test/runtests.jl
index 91aa5d6..3d60b60 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,7 @@
 using StructsOfArrays
 using Test
 using Random
+using StaticArrays
 
 @testset "StructsOfArrays.jl" begin
     @testset "constructor" begin
@@ -52,4 +53,14 @@ using Random
         small[1,1] = TwoField(OneField(1), OneField(2))
         @test small[1,1] == TwoField(OneField(1), OneField(2))
     end
+
+    @testset "StaticArrays" begin
+        type = SArray{Tuple{3},Float64,1,3}
+        regular = rand(MersenneTwister(0), type, 10000)
+        soa = convert(StructOfArrays, regular)
+        @test regular == soa
+        @test sum(regular) ≈ sum(soa)
+        @test eltype(regular) === eltype(soa)
+        @test regular[3] === soa[3]
+    end
 end

From 9c64adcb4ef4440dde69ceeee696f78c8ac8c64e Mon Sep 17 00:00:00 2001
From: Lucas C Wilcox <lucas@swirlee.com>
Date: Wed, 15 Apr 2020 07:29:27 -0700
Subject: [PATCH 4/9] Add a few broadcast tests

---
 test/runtests.jl | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index 3d60b60..a20b9ae 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -39,6 +39,15 @@ using StaticArrays
         @test typeof(similar(small, ComplexF64)) <: StructOfArrays
     end
 
+    @testset "broadcast" begin
+        regular = rand(MersenneTwister(0), ComplexF64, 100)
+        soa = convert(StructOfArrays, regular)
+        @test regular .+ regular == soa .+ soa
+        @test typeof(soa .+ soa) === typeof(soa)
+        @test typeof(regular .+ soa) === typeof(soa)
+        @test sin.(soa)[2] == sin(regular[2])
+    end
+
     @testset "recursive structs" begin
         struct OneField
             x::Int

From a9170061958287dd48d410f73281b64d411142e8 Mon Sep 17 00:00:00 2001
From: Lucas C Wilcox <lucas@swirlee.com>
Date: Tue, 14 Apr 2020 16:01:45 -0700
Subject: [PATCH 5/9] Add basic CuArrays and CUDAnative support

---
 Project.toml           |  6 ++++++
 src/StructsOfArrays.jl | 22 +++++++++++++++++++---
 test/Project.toml      |  3 +++
 test/cuda.jl           | 36 ++++++++++++++++++++++++++++++++++++
 test/runtests.jl       |  7 +++++++
 5 files changed, 71 insertions(+), 3 deletions(-)
 create mode 100644 test/cuda.jl

diff --git a/Project.toml b/Project.toml
index 032ff7c..6ed9601 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,7 +5,13 @@ version = "0.0.3"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+CUDAapi = "3895d2a7-ec45-59b8-82bb-cfc6a382f9b3"
+CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17"
+CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae"
 
 [compat]
 Adapt = "~1"
+CUDAapi = "~4"
+CUDAnative = "~3"
+CuArrays = "~2"
 julia = "≥ 1.3"
diff --git a/src/StructsOfArrays.jl b/src/StructsOfArrays.jl
index 6ccfaf8..bff4341 100644
--- a/src/StructsOfArrays.jl
+++ b/src/StructsOfArrays.jl
@@ -1,6 +1,7 @@
 module StructsOfArrays
 
 export StructOfArrays
+export replace_storage
 
 using Adapt
 
@@ -32,13 +33,28 @@ end
 _type_with_eltype(::Type{<:Array}, T, N) = Array{T, N}
 _type(::Type{<:Array}) = Array
 
-function Adapt.adapt_structure(to, x::StructOfArrays{T, N}) where {T, N}
-    arrays = map(A -> adapt(to, A), x.arrays)
+using CUDAapi
+if has_cuda_gpu()
+    import CuArrays
+    import CuArrays: CuArray
+    _type_with_eltype(::Type{<:CuArray}, T, N) = CuArray{T, N}
+    _type(::Type{<:CuArray}) = CuArray
+
+    import CUDAnative
+    import CUDAnative: CuDeviceArray
+    _type_with_eltype(::Type{<:CuDeviceArray}, T, N) = CuDeviceArray{T, N}
+    _type(::Type{<:CuDeviceArray}) = CuDeviceArray
+end
+
+function replace_storage(f, x::StructOfArrays{T, N}) where {T, N}
+    arrays = map(f, x.arrays)
     TT = typeof(arrays)
-    AT = _type_with_eltype(to, T, N)
+    AT = _type_with_eltype(eltype(arrays), T, N)
     StructOfArrays{T, N, AT, TT}(arrays)
 end
 
+Adapt.adapt_structure(to, x::StructOfArrays{T, N}) where {T, N} = replace_storage(y -> adapt(to, y), x)
+
 function gather_eltypes(T, visited = Set{Type}())
     (!isconcretetype(T) || T.mutable) && throw(ArgumentError("can only create an StructOfArrays of leaf type immutables"))
     if isempty(T.types)
diff --git a/test/Project.toml b/test/Project.toml
index e91b19a..724e558 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,4 +1,7 @@
 [deps]
+CUDAapi = "3895d2a7-ec45-59b8-82bb-cfc6a382f9b3"
+CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17"
+CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StructsOfArrays = "ad108753-388a-4a68-a626-efa2fa9e4278"
diff --git a/test/cuda.jl b/test/cuda.jl
new file mode 100644
index 0000000..bc7be6e
--- /dev/null
+++ b/test/cuda.jl
@@ -0,0 +1,36 @@
+using CuArrays
+using CuArrays: @allowscalar
+using CUDAnative
+CuArrays.allowscalar(false)
+
+@testset "basic" begin
+    type = SArray{Tuple{3},Float64,1,3}
+    N = 1000
+    data = rand(MersenneTwister(0), type, N)
+
+    a = CuArray(data)
+    b = StructOfArrays(type, CuArray, N)
+    c = similar(a)
+    d = replace_storage(CuArray, convert(StructOfArrays, data))
+
+    @test eltype(d) == type
+    @test eltype(eltype(d.arrays)) == Float64
+    @test @allowscalar a[3] === d[3]
+
+    function kernel!(dest, src)
+        i = (blockIdx().x-1)*blockDim().x + threadIdx().x
+        if i <= length(dest)
+            dest[i] = src[i]
+        end
+        return nothing
+    end
+
+    threads = 1024
+    blocks = cld(length(a), threads)
+
+    @cuda threads=threads blocks=blocks kernel!(b, a)
+    @test @allowscalar a[3] === b[3]
+
+    @cuda threads=threads blocks=blocks kernel!(c, b)
+    @test @allowscalar c[3] === b[3]
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index a20b9ae..2f4566d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,6 +2,7 @@ using StructsOfArrays
 using Test
 using Random
 using StaticArrays
+using CUDAapi
 
 @testset "StructsOfArrays.jl" begin
     @testset "constructor" begin
@@ -72,4 +73,10 @@ using StaticArrays
         @test eltype(regular) === eltype(soa)
         @test regular[3] === soa[3]
     end
+
+    if has_cuda_gpu()
+        @testset "CUDA" begin
+            include("cuda.jl")
+        end
+    end
 end

From c4aace83205b1a84d53225e5c3d4fed44fd769fa Mon Sep 17 00:00:00 2001
From: Lucas C Wilcox <lucas@swirlee.com>
Date: Wed, 15 Apr 2020 14:01:19 -0700
Subject: [PATCH 6/9] Add FreeBSD to Travis CI

---
 .travis.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.travis.yml b/.travis.yml
index 29e6ea3..d63dcea 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,6 +1,7 @@
 ## Documentation: http://docs.travis-ci.com/user/languages/julia/
 language: julia
 os:
+  - freebsd
   - linux
   - osx
   - windows

From af009274bc4e8258be5407dd2a7f82541af076e9 Mon Sep 17 00:00:00 2001
From: Lucas C Wilcox <lucas@swirlee.com>
Date: Wed, 15 Apr 2020 17:44:12 -0700
Subject: [PATCH 7/9] Rework README

Here we have updated the README to include the new way of constructing
`StructOfArrays` and have removed dated information.
---
 README.md | 64 ++++++++++++-------------------------------------------
 1 file changed, 14 insertions(+), 50 deletions(-)

diff --git a/README.md b/README.md
index 1f2d721..207d6c4 100644
--- a/README.md
+++ b/README.md
@@ -1,8 +1,5 @@
 # StructsOfArrays
 
-[![Build Status](https://travis-ci.org/simonster/StructsOfArrays.jl.svg?branch=master)](https://travis-ci.org/simonster/StructsOfArrays.jl)
-[![codecov.io](http://codecov.io/github/simonster/StructsOfArrays.jl/coverage.svg?branch=master)](http://codecov.io/github/simonster/StructsOfArrays.jl?branch=master)
-
 A traditional Julia array of immutable objects is an array of structures. Fields
 of a given object are stored adjacent in memory. However, this often inhibits
 SIMD optimizations. StructsOfArrays implements the classic structure of arrays
@@ -13,66 +10,33 @@ contains padding. It is especially useful for arrays of complex numbers.
 
 ## Usage
 
-You can construct a StructOfArrays directly:
+You can construct a `StructOfArrays` directly with:
 
 ```julia
 using StructsOfArrays
-A = StructOfArrays(Complex128, 10, 10)
+A = StructOfArrays(ComplexF64, Array, 10, 10)
 ```
 
-or by converting an AbstractArray:
+or by converting an `AbstractArray`:
 
 ```julia
-A = convert(StructOfArrays, complex(randn(10), randn(10)))
+A = convert(StructOfArrays, rand(ComplexF64, 10, 10))
 ```
 
-Beyond that, there's not much to say. Assignment and indexing works as with
-other AbstractArray types. Indexing a `StructOfArrays{T}` yields an object of
-type `T`, and you can assign objects of type `T` to a given index. The "magic"
-is in the optimizations that the alternative memory layout allows LLVM to
-perform.
-
-While you can create a StructOfArrays of non-`isbits` immutables, this is
-probably slower than an ordinary array, since a new object must be heap
-allocated every time the StructOfArrays is indexed. In practice, StructsOfArrays
-works best with `isbits` immutables such as `Complex{T}`.
-
-## Benchmark
+A `StructOfArrays` can have different storage arrays.  You can construct a
+`CuArray`-based `StructOfArrays` directly with:
 
 ```julia
-using StructsOfArrays
-regular = complex(randn(1000000), randn(1000000))
-soa = convert(StructOfArrays, regular)
-
-function f(x, a)
-    s = zero(eltype(x))
-    @simd for i in 1:length(x)
-        @inbounds s += x[i] * a
-    end
-    s
-end
-
-using Benchmarks
-@benchmark f(regular, 0.5+0.5im)
-@benchmark f(soa, 0.5+0.5im)
+using CuArrays
+A = StructOfArrays(ComplexF64, CuArray, 10, 10)
 ```
 
-The time for `f(regular, 0.5+0.5im)` is:
+or by converting an existing `StructOfArrays` using `replace_storage`:
 
-```
-Average elapsed time: 1.244 ms
-  95% CI for average: [1.183 ms, 1.305 ms]
-Minimum elapsed time: 1.177 ms
-```
-
-and for `f(soa, 0.5+0.5im)`:
-
-```
-Average elapsed time: 832.198 μs
-  95% CI for average: [726.349 μs, 938.048 μs]
-Minimum elapsed time: 713.730 μs
+```julia
+A = replace_storage(CuArray, convert(StructOfArrays, rand(ComplexF64, 10, 10)))
 ```
 
-In this case, StructsOfArrays are about 1.5x faster than ordinary arrays.
-Inspection of generated code demonstrates that `f(soa, a)` uses SIMD
-instructions, while `f(regular, a)` does not.
+Assignment and indexing works as with other `AbstractArray` types. Indexing a
+`StructOfArrays{T}` yields an object of type `T`, and you can assign objects of
+type `T` to a given index.

From 5e49ad5eea2d6401508113f772e6b20498ecf465 Mon Sep 17 00:00:00 2001
From: Lucas C Wilcox <lucas@swirlee.com>
Date: Wed, 15 Apr 2020 16:30:38 -0700
Subject: [PATCH 8/9] Add broadcast support for SoA with CuArray storage

This adds basic support for broadcasting a `StructOfArrays` with `CuArray`
storage.  I could not figure out how to reuse the functionality from
GPUArrays without copying it here.
---
 LICENSE                  |  2 ++
 Project.toml             |  3 +++
 src/StructsOfArrays.jl   | 27 +++++++++-----------
 src/storage/cuda.jl      | 38 ++++++++++++++++++++++++++++
 src/storage/gpuarrays.jl | 53 ++++++++++++++++++++++++++++++++++++++++
 test/Project.toml        |  1 +
 test/cuda.jl             | 23 +++++++++++++++++
 7 files changed, 132 insertions(+), 15 deletions(-)
 create mode 100644 src/storage/cuda.jl
 create mode 100644 src/storage/gpuarrays.jl

diff --git a/LICENSE b/LICENSE
index 29eadaa..9b12cdf 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,5 +1,7 @@
 Copyright (c) 2015: Simon Kornblith.
 Copyright (c) 2018-2019: Valentin Churavy, and other contributors
+Copyright (c) 2016 Simon Danisch
+Copyright (c) 2018 JuliaGPU developers
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
diff --git a/Project.toml b/Project.toml
index 6ed9601..fa87bf1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,10 +8,13 @@ Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 CUDAapi = "3895d2a7-ec45-59b8-82bb-cfc6a382f9b3"
 CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17"
 CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae"
+GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
 Adapt = "~1"
 CUDAapi = "~4"
 CUDAnative = "~3"
 CuArrays = "~2"
+GPUArrays = "~3"
 julia = "≥ 1.3"
diff --git a/src/StructsOfArrays.jl b/src/StructsOfArrays.jl
index bff4341..9b564c1 100644
--- a/src/StructsOfArrays.jl
+++ b/src/StructsOfArrays.jl
@@ -3,7 +3,11 @@ module StructsOfArrays
 export StructOfArrays
 export replace_storage
 
+using Base.Broadcast
+import Base.Broadcast: BroadcastStyle, Broadcasted, AbstractArrayStyle
+
 using Adapt
+using LinearAlgebra
 
 """
     StructOfArrays{T, N, AT, U} <: AbstractArray{T,N}
@@ -33,19 +37,6 @@ end
 _type_with_eltype(::Type{<:Array}, T, N) = Array{T, N}
 _type(::Type{<:Array}) = Array
 
-using CUDAapi
-if has_cuda_gpu()
-    import CuArrays
-    import CuArrays: CuArray
-    _type_with_eltype(::Type{<:CuArray}, T, N) = CuArray{T, N}
-    _type(::Type{<:CuArray}) = CuArray
-
-    import CUDAnative
-    import CUDAnative: CuDeviceArray
-    _type_with_eltype(::Type{<:CuDeviceArray}, T, N) = CuDeviceArray{T, N}
-    _type(::Type{<:CuDeviceArray}) = CuDeviceArray
-end
-
 function replace_storage(f, x::StructOfArrays{T, N}) where {T, N}
     arrays = map(f, x.arrays)
     TT = typeof(arrays)
@@ -162,13 +153,13 @@ end
 end
 
 Base.IndexStyle(::Type{<:StructOfArrays{T,N,A}}) where {T,N,A<:AbstractArray} =  Base.IndexStyle(A)
-Base.BroadcastStyle(::Type{<:StructOfArrays{T,N,A}}) where {T,N,A<:AbstractArray} = Broadcast.ArrayStyle{StructOfArrays{T,N,A}}()
+BroadcastStyle(::Type{<:StructOfArrays{T,N,A}}) where {T,N,A<:AbstractArray} = Broadcast.ArrayStyle{StructOfArrays{T,N,A}}()
 
 function Base.similar(A::StructOfArrays{T1,N,AT}, ::Type{T}, dims::Dims) where {T1,N,AT,T}
     StructOfArrays(T, AT, dims)
 end
 
-function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{StructOfArrays{T1,N,AT}}}, ::Type{T}) where {T1,N,AT,T}
+function Base.similar(bc::Broadcasted{Broadcast.ArrayStyle{StructOfArrays{T1,N,AT}}}, ::Type{T}) where {T1,N,AT,T}
     StructOfArrays(T, AT, Base.to_shape(axes(bc)))
 end
 
@@ -210,4 +201,10 @@ function Base.Array{T,N}(src::StructOfArrays{T,N}) where {T,N}
 end
 Array(src::StructOfArrays{T,N}) where {T,N} = Array{T,N}(src)
 
+include("storage/gpuarrays.jl")
+using CUDAapi
+if has_cuda_gpu()
+    include("storage/cuda.jl")
+end
+
 end
diff --git a/src/storage/cuda.jl b/src/storage/cuda.jl
new file mode 100644
index 0000000..50443f7
--- /dev/null
+++ b/src/storage/cuda.jl
@@ -0,0 +1,38 @@
+import CuArrays
+import CuArrays: CuArray
+_type_with_eltype(::Type{<:CuArray}, T, N) = CuArray{T, N, Nothing}
+_type(::Type{<:CuArray}) = CuArray
+
+import CUDAnative
+import CUDAnative: CuDeviceArray
+_type_with_eltype(::Type{<:CuDeviceArray}, T, N) = CuDeviceArray{T, N}
+_type(::Type{<:CuDeviceArray}) = CuDeviceArray
+
+const AbstractStructOfCuArrays = StructOfArrays{T,N,<:CuArray} where {T, N}
+
+@eval const DestStructOfCuArrays =
+    Union{AbstractStructOfCuArrays,
+          $((:($W where {AT <: AbstractStructOfCuArrays}) for (W, _) in Adapt.wrappers)...),
+          Base.RefValue{<:AbstractStructOfCuArrays} }
+
+function broadcast_kernel!(dest, bc′)
+    i = (CUDAnative.blockIdx().x-1)*CUDAnative.blockDim().x + CUDAnative.threadIdx().x
+    if i <= length(dest)
+        let I = CartesianIndex(CartesianIndices(dest)[i])
+            @inbounds dest[I] = bc′[I]
+        end
+    end
+    return nothing
+end
+
+@inline function Base.copyto!(dest::DestStructOfCuArrays, bc::Broadcasted{Nothing})
+    axes(dest) == axes(bc) || Broadcast.throwdm(axes(dest), axes(bc))
+    bc′ = Broadcast.preprocess(dest, bc)
+
+    threads = 256
+    blocks = cld(length(dest), threads)
+
+    CUDAnative.@cuda threads=threads blocks=blocks broadcast_kernel!(dest, bc′)
+
+    return dest
+end
diff --git a/src/storage/gpuarrays.jl b/src/storage/gpuarrays.jl
new file mode 100644
index 0000000..bb6748e
--- /dev/null
+++ b/src/storage/gpuarrays.jl
@@ -0,0 +1,53 @@
+using GPUArrays
+
+const AbstractStructOfGPUArrays = StructOfArrays{T,N,<:AbstractGPUArray} where {T, N}
+const AbstractStructOfGPUArraysStyle{N} = Broadcast.ArrayStyle{StructOfArrays{T,N,<:AbstractGPUArray}} where T
+
+# Wrapper types otherwise forget that they are StructOfArrays
+for (W, ctor) in Adapt.wrappers
+    @eval begin
+        BroadcastStyle(::Type{<:$W}) where {AT<:StructOfArrays{T,N,A} where {T,N,A}} = BroadcastStyle(AT)
+        backend(::Type{<:$W}) where {AT<:StructOfArrays{T,N,A} where {T,N,A}} = backend(AT)
+    end
+end
+
+# This Union is a hack. Ideally Base would have a
+#     Transpose <: WrappedArray <: AbstractArray
+# and we could define our methods in terms of
+#     Union{AbstractStructOfGPUArrays, WrappedArray{<:Any, <:AbstractStructOfGPUArrays}}
+@eval const DestStructOfGPUArrays =
+    Union{AbstractStructOfGPUArrays,
+          $((:($W where {AT <: AbstractStructOfGPUArrays}) for (W, _) in Adapt.wrappers)...),
+          Base.RefValue{<:AbstractStructOfGPUArrays} }
+
+# Ref is special: it's not a real wrapper, so not part of Adapt,
+# but it is commonly used to bypass broadcasting of an argument
+# so we need to preserve its dimensionless properties.
+BroadcastStyle(::Type{Base.RefValue{AT}}) where {AT<:AbstractStructOfGPUArrays} = typeof(BroadcastStyle(AT))(Val(0))
+backend(::Type{Base.RefValue{AT}}) where {AT<:AbstractStructOfGPUArrays} = backend(AT)
+# but make sure we don't dispatch to the optimized copy method that directly indexes
+function Broadcast.copy(bc::Broadcasted{<:Broadcast.ArrayStyle{StructOfArrays{T,0,<:AbstractGPUArray}}}) where {T}
+    ElType = Broadcast.combine_eltypes(bc.f, bc.args)
+    isbitstype(ElType) || error("Cannot broadcast function returning non-isbits $ElType.")
+    dest = copyto!(similar(bc, ElType), bc)
+    return @allowscalar dest[CartesianIndex()]  # 0D broadcast needs to unwrap results
+end
+
+# Base defines this method as a performance optimization, but we don't know how to do
+# `fill!` in general for all `DestStructOfGPUArrays` so we just go straight to the fallback
+@inline Base.copyto!(dest::DestStructOfGPUArrays, bc::Broadcasted{<:AbstractArrayStyle{0}}) =
+    copyto!(dest, convert(Broadcasted{Nothing}, bc))
+
+## map
+allequal(x) = true
+allequal(x, y, z...) = x == y && allequal(y, z...)
+
+function Base.map!(f, y::DestStructOfGPUArrays, xs::AbstractArray...)
+    @assert allequal(size.((y, xs...))...)
+    return y .= f.(xs...)
+end
+
+function Base.map(f, y::DestStructOfGPUArrays, xs::AbstractArray...)
+    @assert allequal(size.((y, xs...))...)
+    return f.(y, xs...)
+end
diff --git a/test/Project.toml b/test/Project.toml
index 724e558..c97f942 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,6 +2,7 @@
 CUDAapi = "3895d2a7-ec45-59b8-82bb-cfc6a382f9b3"
 CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17"
 CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StructsOfArrays = "ad108753-388a-4a68-a626-efa2fa9e4278"
diff --git a/test/cuda.jl b/test/cuda.jl
index bc7be6e..7304735 100644
--- a/test/cuda.jl
+++ b/test/cuda.jl
@@ -1,6 +1,7 @@
 using CuArrays
 using CuArrays: @allowscalar
 using CUDAnative
+using LinearAlgebra
 CuArrays.allowscalar(false)
 
 @testset "basic" begin
@@ -34,3 +35,25 @@ CuArrays.allowscalar(false)
     @cuda threads=threads blocks=blocks kernel!(c, b)
     @test @allowscalar c[3] === b[3]
 end
+
+@testset "broadcast" begin
+    type = SArray{Tuple{3},Float64,1,3}
+    N = 1000
+    data = rand(MersenneTwister(0), type, N)
+
+    d = replace_storage(CuArray, convert(StructOfArrays, data))
+    e = d .+ d
+    @test @allowscalar e[3] == d[3] .+ d[3]
+    n = norm.(d)
+    @test @allowscalar n[4] == norm(d[4])
+
+    f = convert(StructOfArrays, rand(ComplexF64, 1, 3))
+    g = convert(StructOfArrays, rand(ComplexF64, 3, 1))
+
+    cf = replace_storage(CuArray, f)
+    cg = replace_storage(CuArray, g)
+
+    o = f .+ g
+    co = cf .+ cg
+    @test @allowscalar o[2, 1] == co[2, 1]
+end

From c9a807b292eaf455c23d64594915fba5934b562c Mon Sep 17 00:00:00 2001
From: Lucas C Wilcox <lucas@swirlee.com>
Date: Mon, 20 Apr 2020 17:28:14 -0700
Subject: [PATCH 9/9] Add CUDA example to the README

---
 README.md | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/README.md b/README.md
index 207d6c4..fd1360c 100644
--- a/README.md
+++ b/README.md
@@ -37,6 +37,28 @@ or by converting an existing `StructOfArrays` using `replace_storage`:
 A = replace_storage(CuArray, convert(StructOfArrays, rand(ComplexF64, 10, 10)))
 ```
 
+This array can be used either in a kernel:
+
+```julia
+using CUDAnative
+function kernel!(A)
+    i = (blockIdx().x-1)*blockDim().x + threadIdx().x
+    if i <= length(A)
+        A[i] += A[i]
+    end
+    return nothing
+end
+threads = 256
+blocks = cld(length(A), threads)
+@cuda threads=threads blocks=blocks kernel!(A)
+```
+
+or via broadcasting:
+
+```julia
+A .+= A
+```
+
 Assignment and indexing works as with other `AbstractArray` types. Indexing a
 `StructOfArrays{T}` yields an object of type `T`, and you can assign objects of
 type `T` to a given index.