function as argument of thrust iterator CUDA

Question

I am trying to implement ODEs solver routines running on GPUs using CUDA::Thurst iterators to solve a bunch of equations in the GPU, going to the details, here is a small piece of code:

    #include <thrust/device_vector.h>
    #include <thrust/transform.h> 
    #include <thrust/sequence.h>
    #include <thrust/copy.h> 
    #include <thrust/fill.h>
    #include <thrust/replace.h>
    #include <thrust/functional.h>

    #include <thrust/for_each.h>
    #include <thrust/device_vector.h>
    #include <thrust/iterator/zip_iterator.h>


   #include <iostream>
   #include <math.h>


   __host__ __device__ float f(float x, float y)
   {
     return cos(y)*sin(x);
   }



   struct euler_functor
   {
   const float h;

   euler_functor(float _h) : h(_h) {};

   __host__ __device__
   float operator()( float(*f)(double,double),const float& x, const float& y) const {
   y +=  h * (*f)( x, y );
   x += h;
   }
   };


   int main(void)
   {
   // allocate three device_vectors with 10 elements
   thrust::device_vector<int> X(10);
   // initilaize to random vaues
   thrust::generate(X.begin(), X.end(), rand);
   // apply euler for each element of X
   thrust::for_each(X.begin(),X.end(),euler_functor(f,0.0,X));
   // print the values
   for(int i = 0; i < 10; i++) std::cout<< X[i]<< std::endl;

   }

But when I compile

nvcc euler.cu -o euler.x -lm the following errors occurs:

    lala.cu(29): error: explicit type is missing ("int" assumed)

lala.cu(29): error: expected a ";"

lala.cu(33): error: expression must be a modifiable lvalue

lala.cu(34): error: expression must be a modifiable lvalue

lala.cu(35): warning: missing return statement at end of non-void function "euler_functor::operator()"

lala.cu(46): error: no suitable constructor exists to convert from "float (float, float)" to "euler_functor"

lala.cu(46): error: expected a ")"

it seems like it is not possible use pointers to functions in the way I am trying?

sugestions for better ways to implement the Euler procedure and run it using iterators will be very appreciated.

is the former approach a good compromise between partability and performance?

At the end hopefully the ideal solution for me is be able to define an array of pointer to functions like:

typedef int (*foo_ptr_t)( int );
foo_ptr_t foo_ptr_array[2];

int f1( int );
int f2( int );
foo_ptr_array[0] = f1;
foo_ptr_array[1] = f2;
foo_ptr_array[0]( 1 );

To pass foo_ptr_array as argument to the euler functor. Is it possible?

Thanks for Answer.

Posible improvement:

Is posible define the a set coupled differential equations as fucntors over tuples as I try in following approach? Can I get some error information from the numerical approach to the soution?

It would be

Answer 1

Ultimately, you are asking to take a __device__ function argument in host code, and then pass it as a (function) pointer, in what is ultimately (under the hood) a kernel argument, generated by thrust.

It is illegal to take the address of a __device__ function argument in host code, so passing a __device__ function pointer as an argument this way won't work.

It might be possible to work around this by creating additional __device__ variables (pointers) to store function pointers on the device. Then use cudaGetSymbolAddress to build a table of pointers-to-pointers to functions. This would necessitate running a precursor kernel to set up the function pointers on the device. It seems rather messy.

It might be simpler to parameterize the functor to select a device function based on the parameter. Lke this:

   #include <thrust/device_vector.h>
   #include <thrust/transform.h>
   #include <thrust/sequence.h>
   #include <thrust/copy.h>
   #include <thrust/fill.h>
   #include <thrust/replace.h>
   #include <thrust/functional.h>
   #include <thrust/for_each.h>
   #include <thrust/iterator/zip_iterator.h>

   #include <iostream>
   #include <math.h>


   __host__ __device__ float f1(float x)
   {
     return sinf(x);
   }

   __host__ __device__ float f2(float x)
   {
     return cosf(x);
   }




   struct euler_functor
   {
     unsigned h;

     euler_functor(unsigned _h) : h(_h) {};

     __host__ __device__
     void operator()(float &y) const  {
       if (h == 1) y = f1(y);
       else if (h == 2) y = f2(y);
     }
   };


   int main(void)
   {
     const unsigned N = 8;
     // allocate three device_vectors with 10 elements
     thrust::device_vector<float> X(N);
     // initilaize to random vaues
     thrust::sequence(X.begin(), X.end(),  0.0f, (float)(6.283/(float)N));
     // apply euler for each element of X
     thrust::for_each(X.begin(),X.end(),euler_functor(1));
     // print the values
     for(int i = 0; i < N; i++) std::cout<< X[i]<< std::endl;

     std::cout << "******************" << std::endl;

     thrust::sequence(X.begin(), X.end(),  0.0f, (float)(6.283/(float)N));
     // apply euler for each element of X
     thrust::for_each(X.begin(),X.end(),euler_functor(2));
     // print the values
     for(int i = 0; i < N; i++) std::cout<< X[i]<< std::endl;

   }