Numerical integration by the Romberg method using CUDA -

i tried searching libraries on google numerical integration on cuda couldn't find any.

1) want ask, there libraries available perform integration (of function) on cuda?

2) if write own code on cuda, e.g. implementing romberg integration, how shall proceed? suppose have function, f(x); need calculate integrals of function different intervals e.g. 0.0 - 0.1, ..., 0.2 - 0.3, ..., 1.3 - 2.3? how calculate of them in parallel?

in mind, strategy if have perform, e.g., 1000 integrations, generate 1000 threads, each thread calculates trapzoids error estimates. in case when want calculate trapzoids 1 of integration interval in parallel along other integrals, don't have idea how approach programatically.

as noticed above tera in comment, point of view of parallel programming, integration reduction, simple way implement integration in cuda exploiting primitives of thrust library (see answer simpson's method integrate real valued functions cuda).

below simple example implementing romberg integration method thrust primitives. "direct" translation of corresponding matlab code available @ site, example shows how "simply" matlab codes can ported cuda thurst.

#include <thrust/sequence.h>  #include <thrust/device_vector.h> #include <thrust/host_vector.h>  #define pi_f  3.14159265358979f                 // greek pi in single precision  struct sin_functor {     __host__ __device__     float operator()(float x) const     {         return sin(2.f*pi_f*x);     } };  int main(void) {     int m = 5;                          // --- maximum number of romberg iterations      float     = 0.f;                  // --- lower integration limit     float b     = .5f;                  // --- upper integration limit      float hmin  = (b-a)/pow(2.f,m-1);   // --- minimum integration step size       // --- define matrix romberg approximations , initialize 1.f      thrust::host_vector<float> r(m*m,1.f);      (int k=0; k<m; k++) {          float h = pow(2.f,k-1)*hmin;    // --- step size k-th row of romberg matrix          // --- define integration nodes         int n = (int)((b - a)/h) + 1;         thrust::device_vector<float> d_x(n);         thrust::sequence(d_x.begin(), d_x.end(), a, h);          // --- calculate function values         thrust::device_vector<float> d_y(n);         thrust::transform(d_x.begin(), d_x.end(), d_y.begin(), sin_functor());          // --- calculate integral         r[k*m] = (.5f*h) * (d_y[0] + 2.f*thrust::reduce(d_y.begin() + 1, d_y.begin() + n - 1, 0.0f) + d_y[n-1]);      }      // --- compute k-th column of romberg matrix     (int k=1; k<m; k++) {           // --- matrix of romberg approximations triangular!         (int kk=0; kk<(m-k+1); kk++) {               // --- see romberg integration algorithm             r[kk*m+k] = r[kk*m+k-1] + (r[kk*m+k-1] - r[(kk+1)*m+k-1])/(pow(4.f,k)-1.f);           }       }      // --- define vector rnum numerical approximations     thrust::host_vector<float> rnum(m);      thrust::copy(r.begin(), r.begin() + m, rnum.begin());      (int i=0; i<m; i++) printf("%i %f\n",i,rnum[i]);      getchar();      return 0; }