Can I promote complex number to float number arrays?

Hello! I am new to julia and want to rewrite some of my former C/C++ code in this fascinating language.

I want to solve some equations with complex solution.
As far as I know, most numerical packages only accept real numbers (e.g. NLsolve).

This is not difficult to handle, since we can understand a complex number z as a real array of length 2, so that z=z[1]+im*z[2], and then we can feed this array to the solver.

However, when the dimension of the equations is large enough, one would like to avoid the time of copying the value from a complex array to a real array.
In C/C++, I can do this by a force type conversion, i.e., directly convert a complex pointer to a float pointer, so the corresponding complex array is understood by the compiler as a float array with twice the size.

Is there a way to do the similar thing in julia?

1 Like

You’re looking for reinterpret:

julia> reinterpret(Float64, [1.0+2im, 33])
4-element reinterpret(Float64, ::Vector{ComplexF64}):
  1.0
  2.0
 33.0
  0.0

julia> reinterpret(reshape, Float64, [1.0+2im, 3-4im, 5])
2Ă—3 reinterpret(reshape, Float64, ::Vector{ComplexF64}) with eltype Float64:
 1.0   3.0  5.0
 2.0  -4.0  0.0

If you use ForwardDiff with NLSolve, you will need to do something like reinterpret(complex(eltype(x)), x) rather than specify the exact type.

6 Likes

You CAN’T always reinterpret (“promote”) a number (it has to be an Array of one or more) to an array of two.

If you have a single complex number, it might not be on the heap (and you couldn’t have a “pointer” to it), it might be on the stack yes (and you rather don’t want to have a pointer pointing there), or I guess if could be passed in register[s].

I.e. this works:

julia> x = Complex{Float16}[1.0+2im]
1-element Vector{ComplexF16}:
 Float16(1.0) + Float16(2.0)im

julia> reinterpret(Float16, x)

but:
julia> x = Complex{Float16}(1.0+2im)
Float16(1.0) + Float16(2.0)im

julia> reinterpret(Float16, x)
ERROR: bitcast: expected primitive type value for second argument

Also be careful to match the types:

julia> x = [1.0+2im]; # i.e. ComplexF64

julia> reinterpret(Float16, x)  # to not get nonsensical:
8-element reinterpret(Float16, ::Vector{ComplexF64}):
 0.0
 0.0
 0.0
 1.984
 0.0
 0.0
 0.0
 2.0

A lot of packages are agnostic (if that applies), or should be. I find it intriguing that reinterpreting just worked for you (I realize you can, just wasn’t sure it would always make sense). If/since this does make sense for you (would it always for NLsolve?), I would consider a PR to NLsolve to accept Complex too.

Thank you!
In fact, I have compared the results of NLSolve using complex and equivalent real arrays.
There is no error when using complex array, but the result is very different from the real one.
That’s why I want the reinterprete