Search code examples
c++cudathrust

Using cufft when operating with the thrust library


I want to combine the thrust library and cufft in my project. Thus for testing I wrote

    int length = 5;
    thrust::device_vector<thrust::complex<double> > V1(length);
    thrust::device_vector<cuDoubleComplex> V2(length);
    thrust::device_vector<thrust::complex<double> > V3(length);
    thrust::sequence(V1.begin(), V1.end(), 1);
    thrust::sequence(V2.begin(), V2.end(), 2);
    thrust::transform(V1.begin(), V1.end(), V2.begin(), V3.begin(), thrust::multiplies<thrust::complex<double> >());
    cufftHandle plan;
    cufftPlan1d(&plan, length, thrust::complex<double>, 1);
    cufftExecZ2Z(plan, &V1, &V2, CUFFT_FORWARD);
    for (int i = 0; i < length; i++)
        std::cout << V1[i] << ' ' << V2[i] << ' ' << V3[i] << '\n';
    std::cout << '\n';
    return  EXIT_SUCCESS;

Unfortunately, cufft only accepts arrays such as cuDoubleComplex *a, while thrust::sequence is only working properly with thrust::complex<double>-vectors. When compiling the code above, I get two errors:

error : no operator "=" matches these operands
error : no operator "<<" matches these operands

The first one refers to thrust::sequence(V2.begin(), V2.end(), 2);, while the second one refers to std::cout << V1[i] << ' ' << V2[i] << ' ' << V3[i] << '\n';. If I comment V2, everything works fine.
Is there a conversion between thrust::device_vector<thrust::complex<double>> and cuDoubleComplex *? If not, how can I combine them?


Solution

  • thrust::complex<double> and std::complex<double> share the same data layout with cuDoubleComplex. As a result, all that is required to make your example above work is to cast the data in your device_vector to raw pointers and pass them to cuFFT. Thrust itself can't work with cuDoubleComplex in most operations because that type is a simple container which doesn't define any of the operators which are required to perform any of the operations which Thrust expects for POD types.

    This should work:

    #include <thrust/device_vector.h>
    #include <thrust/transform.h>
    #include <thrust/sequence.h>
    #include <thrust/complex.h>
    #include <iostream>
    #include <cufft.h>
    
    int main()
    {
        int length = 5;
        thrust::device_vector<thrust::complex<double> > V1(length);
        thrust::device_vector<thrust::complex<double> > V2(length);
        thrust::device_vector<thrust::complex<double> > V3(length);
        thrust::sequence(V1.begin(), V1.end(), 1);
        thrust::sequence(V2.begin(), V2.end(), 2);
        thrust::transform(V1.begin(), V1.end(), V2.begin(), V3.begin(), 
                             thrust::multiplies<thrust::complex<double> >());
        cufftHandle plan;
        cufftPlan1d(&plan, length, CUFFT_Z2Z, 1);
        cuDoubleComplex* _V1 = (cuDoubleComplex*)thrust::raw_pointer_cast(V1.data());
        cuDoubleComplex* _V2 = (cuDoubleComplex*)thrust::raw_pointer_cast(V2.data());
    
        cufftExecZ2Z(plan, _V1, _V2, CUFFT_FORWARD);
        for (int i = 0; i < length; i++)
            std::cout << V1[i] << ' ' << V2[i] << ' ' << V3[i] << '\n';
        std::cout << '\n';
        return  EXIT_SUCCESS;
    }