Support for subarrays in linalg

[evelyne-ringoot]

Hello,

I would like to use the triu! and transpose! functions on a non-contiguous view (eg. view(a', 1:2:6,4:2:8)) - is there a way make this possible (ideally for all functions in src/host/linalg.jl; and for copyto! in src/host/abstractarray) without severly increasing runtimes/compiletimes due to multiple-dispatch overhead?

Earlier discussions on this topic:
#452
#458
JuliaGPU/CUDA.jl#1778
JuliaGPU/CUDA.jl#2078

Perhaps some type of a union of subarrays. transposes, and abstractarrays (to avoid switching to AnyGPUArrays; also AnyGPUArrays does not include transposes) ?

Edit: I just saw IndexGPUArray might be an option, if it were expanded with SubArray{T, <:Any, <:LinearAlgebra.Adjoint{T, <:AbstractGPUArray }}

Let me know your thoughts and happy to draft a PR
@maleadt @vchuravy

[maleadt]

As you've noticed before, we're really hesitant to extend everything to AnyGPUArray. Really, we need somebody familiar enough with Julia's array ecosystem to take the time to try and figure this out upstream (e.g., JuliaLang/julia#51910, or some other mechanism to track array storage separate from the array hierarchy)...

In the mean time, a more conservative union that doesn't blow up method complexity (driving up the time to insert them into the method table), increase risk of ambiguities, or include the borderline undefined behavior that Adapt.jl's WrappedArray does, would probably be acceptable too. I guess something sitting in between DenseXArray and StridedXArray, and AnyXArray. Is there any precedent for this in the Julia array ecosystem (I haven't encountered the IndexXArray you're suggesting)?

My only reservation is that the conservativeness of such a union is completely arbitrary; we can index every type of array in a kernel, so the only reason to not have it include every array wrapper is to avoid running in the above issues. So where to you draw the line?

[evelyne-ringoot]

Thanks for the link to the AbstractWrappedArray and appreciate the complexity of that! I found IndexGPUArray in here:

GPUArrays.jl/src/host/indexing.jl

Line 140 in e8e9b03

const IndexGPUArray{T} = Union{AbstractGPUArray{T},

but it might have had a different purpose there?

[maleadt]

Yeah, that's to indicate which arrays can be used as indices, which is unrelated to arrays being indexed within kernels (which is something that's generally possible with any array).

I wonder if at this point we shouldn't consider getting rid of Adapt.jl's horrible WrappedArray altogether and introduce a simpler WrappedGPUArray in GPUArrays only containing the most common Base types, i.e., SubArray, ReshapedArray and ReinterpretArray, doing away with the LinearAlgebra stuff and probably some of the less frequently used wrappers like PermutedDimsArray and LogicalIndex. And if possible the limited two-level wrapping as well, https://github.com/JuliaGPU/Adapt.jl/blob/c80a6ac6c30845fe5fefe694f94b966b5e74af92/src/wrappers.jl#L55-L69.