Support for mathematical operations on complex numbers

[ucanfzh]

When x is a complex-valued CuArray, element-wise operations of basic math function (like exp.(x), sin.(x), etc) gives errors, and Julia crushes for vector operations (sin(x), exp(x), etc).

Original discussion thread on Discourse:
https://discourse.julialang.org/t/support-for-complex-valued-cuarray/16177

[tclements]

Just wondering what would be required to extend CuArray to support transcendental operations with complex type?

[maleadt]

Which functions do you need specifically? Basically, we now rely on libdevice, a bitcode library provided by NVIDIA, for most mathematical operations. None of these support complex numbers. We'd need to carefully check which existing Base functions are GPU compatible nowadays, and/or provide alternatives. I did a quick check and neither sin nor exp, both of which are implemented in pure Julia, seem to work out of the box:

julia> using CuArrays

julia> a = CuArrays{Complex{Int}}(undef, (2,2))
ERROR: TypeError: in Type{...} expression, expected UnionAll, got Module
Stacktrace:
 [1] top-level scope at none:0

julia> a = CuArray{Complex{Int}}(undef, (2,2))
2×2 CuArray{Complex{Int64},2}:
 0+0im  0+0im
 0+0im  0+0im

# basic stuff just works

julia> a .+ 1
2×2 CuArray{Complex{Int64},2}:
 1+0im  1+0im
 1+0im  1+0im

# some functions don't

julia> sin.(a)
ERROR: MethodError: no method matching sin(::Complex{Int64})
You may have intended to import Base.sin

julia> using CUDAnative

julia> CUDAnative.sin(x::Complex) = Base.sin(x)

ERROR: InvalidIRError: compiling #23(CuArrays.CuKernelState, CuDeviceArray{Complex{Float64},2,CUDAnative.AS.Global}, Base.Broadcast.Broadcasted{Nothing,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}},typeof(CUDAnative.sin),Tuple{Base.Broadcast.Extruded{CuDeviceArray{Complex{Int64},2,CUDAnative.AS.Global},Tuple{Bool,Bool},Tuple{Int64,Int64}}}}) resulted in invalid LLVM IR

# nope

julia> CUDAnative.exp(x::Complex) = Base.exp(x)

julia> exp.(a)
ERROR: LLVM error: Cannot select: 0x52eacf0: i64,glue = sube Constant:i64<0>, 0x52eac20, 0x52eac88:1

# nope

[tclements]

It would be nice to have sin, cos and exp at least. Most other functions (sinh, cosh, etc) are combinations of sin, cos and exp. My use case is spectral analysis, which is wholly dependent on manipulating complex numbers. Here is my hackish workaround for now:

julia> using GPUArrays

julia> using CuArrays

julia> import Base.sin

julia> function sin(A::GPUArray{ComplexF64})
             return sin.(real(A)) .* cosh.(imag(A))  .+ im .* cos.(real(A)) .* sinh.(imag(A))
       end
sin (generic function with 15 methods)

julia> a = cu(rand(ComplexF64,2))
2-element CuArray{Complex{Float64},1}:
 0.6257075082417543 + 0.9629084552943576im
 0.4461085482428031 + 0.6315414590805859im

julia> sin(a)
2-element CuArray{Complex{Float64},1}:
 0.8788238879768442 + 0.9068108257247335im
 0.5203986812134053 + 0.6083694003870889im

julia> isapprox(collect(sin(a)), sin.(collect(a)))
true

[ucanfzh]

It would be nice also to support fft,ifft, etc for Fourier analysis. Thanks.

…

On Fri, 15 Feb 2019 at 00:56, Tim Clements ***@***.***> wrote: It would be nice to have sin, cos and exp at least. Most other functions ( sinh, cosh, etc) are combinations of sin, cos and exp. My use case is spectral analysis, which is wholly dependent on manipulating complex numbers. Here is my hackish workaround for now: julia> using GPUArrays julia> using CuArrays julia> import Base.sin julia> function sin(A::GPUArray{ComplexF64}) return sin.(real(A)) .* cosh.(imag(A)) .+ im .* cos.(real(A)) .* sinh.(imag(A)) end sin (generic function with 15 methods) julia> a = cu(rand(ComplexF64,2)) 2-element CuArray{Complex{Float64},1}: 0.6257075082417543 + 0.9629084552943576im 0.4461085482428031 + 0.6315414590805859im julia> sin(a) 2-element CuArray{Complex{Float64},1}: 0.8788238879768442 + 0.9068108257247335im 0.5203986812134053 + 0.6083694003870889im julia> isapprox(collect(sin(a)), sin.(collect(a))) true — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#165 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AVtheEL9O2XwIIONwi_KXZqo35Ix8P2Hks5vNZUtgaJpZM4XbKbB> .

[coezmaden]

Have there been any updates on this issue? Would be helpful in my research as well.

As of 14.04.2020 exp seems to work fine.

julia> a = cu(rand(Complex{Float32}, 2))

2-element CuArray{Complex{Float32},1,Nothing}:
 0.48850453f0 + 0.341748f0im
 0.39827824f0 + 0.9656277f0im

# declaring a function for element-wise exponentiatation of a C^n vector
julia> function expvector!(a::CuArray{ComplexF32})
               a = exp.(a)
               return a
       end

expvector! (generic function with 1 method)

julia> expvector!(a)
2-element CuArray{Complex{Float32},1,Nothing}:
  1.5356216f0 + 0.546228f0im
 0.84724027f0 + 1.2247753f0im

#check wolfram : true

An implementation of sin and cos for complex numbers exists in Julia Base https://github.com/JuliaLang/julia/blob/master/base/complex.jl#L804
However, I couldn´t get it to run.

julia> function sinvector!(a::CuArray{ComplexF32})
               a = Base.sin.(a)
               return a
       end

ERROR: MethodError: no method matching sin(::Complex{Float32})
You may have intended to import Base.sin

[maleadt]

Broadcasting automatically changes invocations of Base.sin to CUDAnative.sin. PRs for extended type coverage are welcome.

[maleadt]

Should be fixed in the back-end, e.g., JuliaGPU/CUDA.jl#750.