Question about how to use on gpu

[yolhan83]

Hello, I wonder why this is not working on gpu :

julia> using CUDA,FiniteDiff
julia> f(x) = sum(abs2,x)
f (generic function with 1 method)

julia> x = rand(100,100) |> CuArray;

julia> f(x)
3298.8592709432687

julia> FiniteDiff.finite_difference_gradient(f,x);
ERROR: MethodError: no method matching compute_epsilon(::Val{:central}, ::CuArray{Float64, 2, CUDA.DeviceMemory}, ::Float64, ::Float64, ::Bool)
The function `compute_epsilon` exists, but no method is defined for this combination of argument types.

julia> FiniteDiff.finite_difference_gradient(f,collect(x)); # works

[ChrisRackauckas]

I don't currently have a GPU to debug this: can you share the full stack trace so I can see if I can eyeball it?

[yolhan83]

Yes, here it is

MethodError: no method matching compute_epsilon(::Val{:central}, ::CuArray{Float64, 2, CUDA.DeviceMemory}, ::Float64, ::Float64, ::Bool)
The function `compute_epsilon` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  compute_epsilon(::Val{:central}, ::T, ::Real, ::Real, ::Any) where T<:Number
   @ FiniteDiff C:\Users\yolha\.julia\packages\FiniteDiff\vund7\src\epsilons.jl:13
  compute_epsilon(::Val{:hcentral}, ::T, ::Real, ::Real, ::Any) where T<:Number
   @ FiniteDiff C:\Users\yolha\.julia\packages\FiniteDiff\vund7\src\epsilons.jl:17
  compute_epsilon(::Val{:central}, ::T, ::Real, ::Real) where T<:Number
   @ FiniteDiff C:\Users\yolha\.julia\packages\FiniteDiff\vund7\src\epsilons.jl:13
  ...


Stacktrace:
 [1] #finite_difference_gradient!#9
   @ C:\Users\yolha\.julia\packages\FiniteDiff\vund7\src\gradients.jl:280 [inlined]
 [2] finite_difference_gradient!
   @ C:\Users\yolha\.julia\packages\FiniteDiff\vund7\src\gradients.jl:228 [inlined]
 [3] finite_difference_gradient(f::typeof(f), x::CuArray{Float64, 2, CUDA.DeviceMemory}, fdtype::Val{:central}, returntype::Type, inplace::Val{true}, fx::Nothing, c1::Nothing, c2::Nothing; relstep::Float64, absstep::Float64, dir::Bool)
   @ FiniteDiff C:\Users\yolha\.julia\packages\FiniteDiff\vund7\src\gradients.jl:163
 [4] finite_difference_gradient(f::Function, x::CuArray{Float64, 2, CUDA.DeviceMemory}, fdtype::Val{:central}, returntype::Type, inplace::Val{true}, fx::Nothing, c1::Nothing, c2::Nothing)
   @ FiniteDiff C:\Users\yolha\.julia\packages\FiniteDiff\vund7\src\gradients.jl:129
 [5] top-level scope
   @ In[12]:3

meaning, it enters here in the code,

function finite_difference_gradient!(
    df,
    f,
    x,
    cache::GradientCache{T1,T2,T3,T4,fdtype,returntype,inplace};
    relstep=default_relstep(fdtype, eltype(x)),
    absstep=relstep,
    dir=true) where {T1,T2,T3,T4,fdtype,returntype,inplace}

    # NOTE: in this case epsilon is a vector, we need two arrays for epsilon and x1
    # c1 denotes x1, c2 is epsilon
    fx, c1, c2, c3 = cache.fx, cache.c1, cache.c2, cache.c3
    if fdtype != Val(:complex) && ArrayInterface.fast_scalar_indexing(c2)
        @. c2 = compute_epsilon(fdtype, x, relstep, absstep, dir)
        copyto!(c1, x)
    end
    copyto!(c3, x)
    if fdtype == Val(:forward)
        ......
    elseif fdtype == Val(:central)
        @inbounds for i ∈ eachindex(x)
            if ArrayInterface.fast_scalar_indexing(c2)
                epsilon = ArrayInterface.allowed_getindex(c2, i) * dir
            else
                epsilon = compute_epsilon(fdtype, x, relstep, absstep, dir) * dir # !!!!!!!!!!!!!!!!!!!!!!!!!!!
            end

and breaks because x should be a number and its the array that enters.
Also, it's the same for :forward

[ChrisRackauckas]

A new patch should come out in an hour. That should handle this. compute_epsilon(fdtype, x, relstep, absstep, dir) should've been compute_epsilon(fdtype, one(eltype(x)), relstep, absstep, dir). I don't have a good way to test it, because you really only hit this with GPU arrays, but it's at least a clear typo so I put the patch out for now.

[yolhan83]

Thank you, the gradient pass now but is completly false,

 using CUDA,FiniteDiff
f(x) = sum(abs2,x)
x = (Float32[1 2;3 4]) |> CuArray;
display(FiniteDiff.finite_difference_gradient(f,x))
FiniteDiff.finite_difference_gradient(f,collect(x))

results :

2×2 CuArray{Float32, 2, CUDA.DeviceMemory}:
 -3046.81  -3045.81
 -3044.81  -3043.81
2×2 Matrix{Float32}:
 2.00014  3.99989
 6.00023  7.99998

and for forward,

 using CUDA,FiniteDiff
f(x) = sum(abs2,x)/length(x)
x = (Float32[1 2;3 4]) |> CuArray;
display(FiniteDiff.finite_difference_gradient(f,x,Val(:forward)))
FiniteDiff.finite_difference_gradient(f,collect(x),Val(:forward))

results :

2×2 CuArray{Float32, 2, CUDA.DeviceMemory}:
 -86889.3  -86889.3
 -86889.3  -86889.3
2×2 Matrix{Float32}:
 1.99979  3.99957
 5.99936  7.99915

It works fine when x is a vector but not when it is an array ex :

display(FiniteDiff.finite_difference_gradient(f,x[:]))
FiniteDiff.finite_difference_gradient(f,collect(x[:]))

results :

4-element CuArray{Float32, 1, CUDA.DeviceMemory}:
 2.000141
 6.000035
 3.999991
 7.999982
4-element Vector{Float32}:
 2.000141
 6.000035
 3.999991
 7.999982

[ChrisRackauckas]

@oscardssmith can you look into this? I don't have a GPU and it needs to look at the no-indexing GPU branches.