State of GPU support

I wanted to ask for an overview of the current state of GPU support of this package. It appears as though there are several issues related to whether the package works nicely on the GPU (i.e. with GPU arrays as inputs) and several proposed solutions, but getting KernelFunctions.jl to work on the GPU seems to be delayed by other things breaking, like AD.

I was wondering whether a clear path forward is already emerging. Since I've done some GPU work before, I'd be happy to help getting this package work on the GPU.

For what I mean by GPU support, at least the following should be possible:

using CUDA, KernelFunctions
CUDA.allowscalar(false)
x = CUDA.rand(16)
k = SEKernel()
kernelmatrix(k, x)

Start by reviewing the related GPU-array and automatic-differentiation issues referenced in the discussion, then run the shown Julia example with CUDA, KernelFunctions, SEKernel, kernelmatrix, and CUDA.allowscalar(false). Done means the example works with a CUDA array input and a clear path forward for GPU support is documented.