Complex number support in OpenCL

Go To StackoverFlow.com

4

I know that OpenCL doesn't have support for complex numbers, and by what I've read this feature isn't going to show up anytime soon. Still, several examples make use of complex numbers in OpenCL kernels (FFT implementations for instance).

Does anybody have experience with this? What would be the "best" method to enable complex number support in OpenCL? I'd assume using a float2 to contain the real and imaginary parts, but should I write a set of macros, or are inline functions better? Does anybody know if a set of functions/macros already exists for this purpose?

2012-04-04 17:16
by Emanuel Ey
The article here illustrates how to do it: http://developer.amd.com/resources/documentation-articles/articles-whitepapers/opencl-optimization-case-study-fast-fourier-transform-part-1 - ChrisF 2014-10-01 16:15


6

So, as I needed a set of functions to handle complex numbers in OpenCL I ended up implementing a set of them. Specifically i needed sum and subtraction (trivial, can be done with standard vector operations), multiplication, division, getting a complex's modulus, argument (or angle) and square root.
relevant wikipedia articles:
http://en.wikipedia.org/wiki/Complex_number#Absolute_value_and_argument
http://en.wikipedia.org/wiki/Square_root#Principal_square_root_of_a_complex_number
This is mostly trivial, but it does take some time, so in the hope in might save someone this time, here goes:

//2 component vector to hold the real and imaginary parts of a complex number:
typedef float2 cfloat;

#define I ((cfloat)(0.0, 1.0))


/*
 * Return Real (Imaginary) component of complex number:
 */
inline float  real(cfloat a){
     return a.x;
}
inline float  imag(cfloat a){
     return a.y;
}

/*
 * Get the modulus of a complex number (its length):
 */
inline float cmod(cfloat a){
    return (sqrt(a.x*a.x + a.y*a.y));
}

/*
 * Get the argument of a complex number (its angle):
 * http://en.wikipedia.org/wiki/Complex_number#Absolute_value_and_argument
 */
inline float carg(cfloat a){
    if(a.x > 0){
        return atan(a.y / a.x);

    }else if(a.x < 0 && a.y >= 0){
        return atan(a.y / a.x) + M_PI;

    }else if(a.x < 0 && a.y < 0){
        return atan(a.y / a.x) - M_PI;

    }else if(a.x == 0 && a.y > 0){
        return M_PI/2;

    }else if(a.x == 0 && a.y < 0){
        return -M_PI/2;

    }else{
        return 0;
    }
}

/*
 * Multiply two complex numbers:
 *
 *  a = (aReal + I*aImag)
 *  b = (bReal + I*bImag)
 *  a * b = (aReal + I*aImag) * (bReal + I*bImag)
 *        = aReal*bReal +I*aReal*bImag +I*aImag*bReal +I^2*aImag*bImag
 *        = (aReal*bReal - aImag*bImag) + I*(aReal*bImag + aImag*bReal)
 */
inline cfloat  cmult(cfloat a, cfloat b){
    return (cfloat)( a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x);
}


/*
 * Divide two complex numbers:
 *
 *  aReal + I*aImag     (aReal + I*aImag) * (bReal - I*bImag)
 * ----------------- = ---------------------------------------
 *  bReal + I*bImag     (bReal + I*bImag) * (bReal - I*bImag)
 * 
 *        aReal*bReal - I*aReal*bImag + I*aImag*bReal - I^2*aImag*bImag
 *     = ---------------------------------------------------------------
 *            bReal^2 - I*bReal*bImag + I*bImag*bReal  -I^2*bImag^2
 * 
 *        aReal*bReal + aImag*bImag         aImag*bReal - Real*bImag 
 *     = ---------------------------- + I* --------------------------
 *            bReal^2 + bImag^2                bReal^2 + bImag^2
 * 
 */
inline cfloat cdiv(cfloat a, cfloat b){
    return (cfloat)((a.x*b.x + a.y*b.y)/(b.x*b.x + b.y*b.y), (a.y*b.x - a.x*b.y)/(b.x*b.x + b.y*b.y));
}


/*
 *  Square root of complex number.
 *  Although a complex number has two square roots, numerically we will
 *  only determine one of them -the principal square root, see wikipedia
 *  for more info: 
 *  http://en.wikipedia.org/wiki/Square_root#Principal_square_root_of_a_complex_number
 */
 inline cfloat csqrt(cfloat a){
     return (cfloat)( sqrt(cmod(a)) * cos(carg(a)/2),  sqrt(cmod(a)) * sin(carg(a)/2));
 }
2012-04-06 10:45
by Emanuel Ey
That one-letter #define I is going to bite you hard one day ;-) - rubenvb 2014-06-26 08:54
@rubenvb why is that? do you have an example - Emanuel Ey 2014-07-21 12:50
You can substitute the carg(cfloat a) with the atan2 function that is already implemented in OpenCL: https://www.khronos.org/registry/cl/sdk/1.0/docs/man/xhtml/atan.htm - Dirklinux 2015-10-20 08:57


4

PyOpenCL has a somewhat more complete and robust implementation of complex numbers in OpenCL:

https://github.com/pyopencl/pyopencl/blob/master/pyopencl/cl/pyopencl-complex.h

2012-04-18 17:38
by Andreas Klöckner
Guter Tip, danke :) From your repository "Note that addition (real + complex) and multiplication (complexcomplex) are defined, but yield wrong results". what do you mean by this? real+complex addition was actually not defined and complexcomplex multiplication will expand to return (TP)( a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x); } which looks fine to me or is there something I'm overlooking - Emanuel Ey 2012-04-19 14:07
link update: https://github.com/pyopencl/pyopencl/blob/master/pyopencl/cl/pyopencl-complex. - tesch1 2015-03-18 16:59
Ads