thrust stuff silently fails when using thrust::backend::vector instead of thrust::device_vector

[slayoo]

Here's an example code:

#define THRUST_DEVICE_SYSTEM THRUST_DEVICE_SYSTEM_CPP

#include <thrust/system/cpp/vector.h>

#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/external/thrust/thrust_algebra.hpp>
#include <boost/numeric/odeint/external/thrust/thrust_operations.hpp>
#include <boost/numeric/odeint/external/thrust/thrust_resize.hpp>

struct rhs 
{
  void operator()(
    const thrust::cpp::vector<float> &psi,
    thrust::cpp::vector<float> &dot_psi,
    const float /* t */
  )  
  {
    std::cerr << "size: " << dot_psi.size() << std::endl;
  }  
};

int main() 
{ 
  boost::numeric::odeint::euler<
    thrust::cpp::vector<float>, // state_type
    float,                      // value_type
    thrust::cpp::vector<float>, // deriv_type
    float,                      // time_type
    boost::numeric::odeint::thrust_algebra,
    boost::numeric::odeint::thrust_operations 
  > chem_stepper;

  thrust::cpp::vector<float> v(11);
  chem_stepper.do_step(rhs(), v, 0, 1);
}

Compilation goes fine, but when one tries to use it, there is a segfault due to internal stepper variables not being properly resized:

$ g++ bug.cpp 
$ ./a.out 
size: 0
Segmentation fault

Using cpp::vector is convenient (if not the only possible way) when using multiple Thrust backends in a single translation unit.

A solution would be to repeat all the definitions in boost/numeric/odeint/external/thrust/*.hpp not only for host_vector and device_vector, but also for thrust::cpp::vector, thrust::cuda::vector, thrust::tbb::vector and thrust::omp::vector.

HTH, Sylwester

[mariomulansky]

All of those vectors can be used in the thrust::for_each algorithms i suppose? So we would use the same thrust_algebra for all of them?

[slayoo]

Yes. Thanks for looking into it!

[ddemidov]

It seems that all kinds of vectors in thrust inherit from thrust::detail::vector_base<T,A>. This could be used to create generic thrust vector specifications for odeint.

[mariomulansky]

good idea denis, this would make things rather simple - i will look into this

[mariomulansky]

I figured an easy approach is to write some macros. I can't test the thrust::cpp::vector because I only have cuda 5 installed where there is no cpp::vector it seems. However - maybe you can give it a try and add the macro calls to external/thrust/thrust_resize.hpp ?

[slayoo]

Thanks Mario!

I confirm it works with cpp, I'll test tomorrow with CUDA and OpenMP

There is probably a need to enclose the ::vector stuff in a proprocessor condition on Thrust version. Here's a relevant Thrust commit that intorduced the ::vector files:

https://github.com/thrust/thrust/commit/9a77b6542770071eb67aefde68ea90d531ac1253

but I haven't yet figured out in which release it appeared.

Also the thrust_resize.hpp will have to include something like:

#include <thrust/system/cpp/vector.h>
#include <thrust/system/omp/vector.h>
#include <thrust/system/tbb/vector.h>
#if defined(__NVCC__)
#  include <thrust/system/cuda/vector.h> 
#endif

... also only if thrust version is enough.

Thanks, Sylwester

[mariomulansky]

I'm leaving into holiday tomorrow and I will not be able to put more work in it before September, sorry. But I hope for the time being the Macros are making it easy enough for you to continue your work?

[slayoo]

Yes, it definitely helps. Thanks! I've added a pull request with the code that uses the new macros . Tested with cpp and omp backends using the code I'm workin on, not tested with odeint examples.

HTH, Sylwester

[mariomulansky]

Thanks a lot slayoo, I've merged the pull request and added the corresponding algebra_dispatcher and operations_dispatcher specializations. Maybe you can have a look if your stuff still works nicely. Now you should be able to use the thrust backend vectors also without explicitely defining the thrust_algebra or thrust_operations, i.e.

boost::numeric::odeint::euler<
    thrust::cpp::vector<float>, // state_type
    float,                      // value_type
  > chem_stepper;

[mariomulansky]

was the fix sufficient, can this be closed?

[slayoo]

Yes. I've just double checked with the current snapshot of odeint, the above short example as well as with our code that caused the problem originally. Thanks!

What about adding the above example as a unit test? (sorry, I haven't checked if it is not already done)

Sylwester