Slow complex exponential evaluation

Specific Domains Numerics
1

I am having some issues with code and I can’t seem to reach a solution. Specifically I am re-writing some Matlab code and for one or two things I can’t seem to match the speed of Matlab.

Matlab Code:

a=[0:2^24-1]
b=exp(2*pi*a*i)

This is pretty fast but I can’t seem to match the execution in Julia:

a=collect(0:24-1)
b=exp(2*pi*im*a)

This is very, very slow ~10s after compilation. I tried to make an explicit loop and attempt that and only saw negligible improvement:

function imexp_1d(x::Array)

        out::Array{Complex{Float64}}=Array(Complex{Float64},size(x))
        for i in eachindex(x)
                out[i]=cos(x[i]*2.0*pi)+sin(x[i]*2.0*pi)*im
                #out[i]=complex(cos(x[i]*2.0*pi),sin(x[i]*2.0*pi))
        end
        return out
end
tt=collect(Float64,0:2^24-1)
@time imexp_1d(tt);
@time imexp_1d(tt);

The output is still ~9s. Am I missing something? Any help would be appreciated.

2

that’s not strictly typed. An Array has two paramters: the element type and the dimension. Instead, do

out=similar(x,Complex{Float64})

There’s no reason to declare the type either.

Edit

Nevermind, I never knew you could do that with the Array constructor. Huh, that’s cool. I still use similar.

3

This is from matlab executing the operation with multiple threads.

4

For reference, here’s how I would do it:

function imexp{T}(x::AbstractArray{T})
  out = similar(x,Complex{T})
  imexp!(out,x)
  out
end
function imexp!(out,x::AbstractArray)
  for i in eachindex(x)
    out[i]=cos(x[i]*2pi)+sin(x[i]*2pi)*im
  end
end
tt=collect(Float64,0:2^24-1)
@time imexp(tt);
@time imexp(tt);

I get:

5.993910 seconds (6 allocations: 256.000 MB, 1.38% gc time)
5.971867 seconds (6 allocations: 256.000 MB, 1.15% gc time)

If you want threads, just do

function imexp!(out,x::AbstractArray)
  Theads.@threads for i in eachindex(x)
    out[i]=cos(x[i]*2.0*pi)+sin(x[i]*2.0*pi)*im
  end
end

and make sure you enabled multithreading:

https://docs.julialang.org/en/stable/manual/parallel-computing/#setup

Edit

I realized that the Array constructor usage you did there actually works, so there’s no real timing difference. So it must all be multithreading.

1 Like
5

Also I assume you mean 2^24? It’ll be very surprising if doing exp(::Complex128) on 24 elements takes 10s for any machine we support…

6

cis

1 Like
7

Thanks! Yea i guess it must be the multithreading. I just tried your modified code with 4 threads going and i got ~2.5s timing which is much better. I thought it might be the multithreading thing but when I looked it up Matlab’s docs said that they only used multhreading automatically for LinAlg type operations.

Also, is there a particular reason for using the imexp!() separately? Does it have to do with type-stability?

8

That’s an in-place operation. If your output array is already allocated, it can be quicker depending on the size of the array, and also it can put less of a load in the garbage collector if you’re using this in a loop.

1 Like
9

Gotcha! Thanks a bunch.