support for complex values about cusp-library HOT 12 CLOSED

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

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.

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

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.

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

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.

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()).

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;
}

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

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

This issue was closed by revision 8e7d10a16b.

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!