Basic math operations with nested adjoint and transpose

[albertomercurio]

Hello,

I need to make the following code working for CUDA too:

A = rand(ComplexF64, 10, 10)
adjoint(transpose(A')) - adjoint(transpose(A'))

10×10 Matrix{ComplexF64}:
 0.0+0.0im  0.0+0.0im  0.0+0.0im  …  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  …  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im

But if I try to do it with a CuArray if fails:

A = CUDA.rand(ComplexF64, 10, 10)
adjoint(transpose(A')) - adjoint(transpose(A'))

Scalar indexing is disallowed.
Invocation of getindex resulted in scalar indexing of a GPU array.
This is typically caused by calling an iterating implementation of a method.
Such implementations *do not* execute on the GPU, but very slowly on the CPU,
and therefore are only permitted from the REPL for prototyping purposes.
If you did intend to index this array, annotate the caller with @allowscalar.

It would be nice to introduce such an implementation, perhaps even for CuSparse matrices.

[maleadt]

That's unfortunately not straightforward. Closing as a dup of JuliaGPU/Adapt.jl#21.

For some more context, this needs better support in Base, either through a better type-hierarchy or via a trait-like system. We cannot simply add specializations for every combination of nested wrappers here in CUDA.jl, not only would that take many definitions, it would also severely impact load time.