support for complex values

[GoogleCodeExporter]

following the implementation of std::complex<T>

Specializations for float and double should use float2 and double2 for 
coalescing.

Original issue reported on code.google.com by wnbell on 31 Jul 2010 at 6:12

[GoogleCodeExporter]

I think it's important to make the complex type binary compatible with 
cuComplex/cuDoubleComplex.

Original comment by filipe.c...@gmail.com on 29 Sep 2010 at 7:51

[GoogleCodeExporter]

Attached is my first proposal for complex numbers.
I didn't implement long doubles, nor the >> operator.
I have only tested it a little bit at the moment, but I'm planning to test it 
more.

Original comment by filipe.c...@gmail.com on 1 Oct 2010 at 1:37

Attachments:

• cusp_complex.h

[GoogleCodeExporter]

I wonder if it would not be better to add complex numbers in thrust instead of 
cusp.

Original comment by filipe.c...@gmail.com on 1 Oct 2010 at 5:53

[GoogleCodeExporter]

Thanks Fillipe, this is a valuable contribution!

I think we will eventually add complex numbers to Thrust, just not immediately. 
 Our plan is to create a thrust::numerical:: sub-library that will provide 
functionality for numerical computations.  For example, it will provide a 
container like cusp::array1d that does not initialize its contents, common 
vector types like float2 and int4, and the standard cmath functions among other 
things.

Anyway, we would be happy to incorporate complex into Cusp as soon as you're 
done.  And once thrust::numerical:: exists we'll move the implementation over 
there.

Original comment by wnbell on 1 Oct 2010 at 6:08

[GoogleCodeExporter]

One more thing :)

We'll want to make cusp::complex (and eventually thrust::numerical::complex) 
conform to the definition in the C++ standard [1] (see Section 26.4 on page 
886).  I think your implementation is on track, it just needs a few minor 
changes (e.g. needs to define complex<T>::value_type).

[1] http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2010/n3092.pdf

Original comment by wnbell on 1 Oct 2010 at 6:22

[GoogleCodeExporter]

I've reached a point where I think cusp::complex is usable.
I've attached here the code (also available from my cloned repository 
(https://code.google.com/r/filipecmaia-cusp-complex/).

Any feedback appreciated.

Original comment by filipe.c...@gmail.com on 5 Oct 2010 at 10:54

Attachments:

• complex.h
• complex.cu

[GoogleCodeExporter]

Doing some more tests made me realize of the need to overload the non member 
functions for the non complex types otherwise code in the cusp namespace that 
tries to call for example sqrt() will get the cusp::sqrt() from complex.h which 
does not accept non complex arguments.
To fix this problem I created a cmath.h file which simply redirects all such 
calls to the appropriate std:: function (for example std::sqrt()).

Original comment by filipe.c...@gmail.com on 7 Oct 2010 at 11:24

Attachments:

• cmath.h

[GoogleCodeExporter]

Thanks for your continued work on this Filipe! I'll pull the changes from your 
HG clone and merge them with mainline Cusp within the next few days.

One suggestion:

Instead of using __CUDA_ARCH__ to determine whether to test doubles you should 
detect that at runtime with something like the following code.  The reason this 
approach is better is that __CUDA_ARCH__ only makes sense if the user specifies 
exactly one -arch during compilation.  If the user specifies multiple arch 
targets with the -gencode option they'll get undefined behavior.

bool device_supports_double(void)
{
    int current_device = -1;
    cudaDeviceProp properties;

    cudaError_t error = cudaGetDevice(¤t_device);
    if(error)
        throw thrust::system_error(error, thrust::cuda_category());

    if(current_device < 0)
        throw thrust::system_error(cudaErrorNoDevice, thrust::cuda_category());

    // the properties weren't found, ask the runtime to generate them
    error = cudaGetDeviceProperties(&properties, current_device);

    if(error)
      throw thrust::system_error(error, thrust::cuda_category());

    return properties.major >= 1 && properties.minor >= 3;
}

Original comment by wnbell on 7 Oct 2010 at 11:59

[GoogleCodeExporter]

In short, using __CUDA_ARCH__ in host code is potentially dangerous.  Using it 
within device code is perfectly fine however.

Original comment by wnbell on 8 Oct 2010 at 12:01

[GoogleCodeExporter]

Thanks for the tip. It's (hopefully) fixed now.

Original comment by filipe.c...@gmail.com on 8 Oct 2010 at 12:35

[GoogleCodeExporter]

This issue was closed by revision 8e7d10a16b.

Original comment by wnbell on 17 Oct 2010 at 8:45

• Changed state: Fixed

[GoogleCodeExporter]

Sorry for the delay Filipe, I just got around to merging your complex changes 
into the trunk.  Since then I made a handful of changes so cusp::complex 
inter-operates better with scalar values and std::complex.  I plan to refactor 
complex.h a little to remove some of the redundancy between complex<float> and 
complex<double>.  After that I'll need to add some additional tests to make 
sure the various matrix operations support complex values, but with any luck 
those should "just work".  Once we're comfortable that everything works as 
expected I'll make an announcement on cusp-users.

Thanks again!

Original comment by wnbell on 17 Oct 2010 at 8:57