Romberg's integration of real valued functions by CUDA Thrust

From the point of view of parallel programming, integration is basically a reduction, so that a very simple way to implement integration in CUDA is exploiting the primitives of the Thrust library.

Below is a simple example implementing the Romberg integration method by the Thrust primitives.
It is a “direct” translation of the corresponding Matlab code available at this site, so this example also shows how “simply” some Matlab codes can be ported to CUDA by 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 a     = 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 the matrix for Romberg approximations and initialize to 1.f

   thrust::host_vector<float> R(M*M,1.f);

   for (int k=0; k<M; k++)
      float h = pow(2.f,k-1)*hmin;      // --- Step size for the k-th row of the 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 the k-th column of the Romberg matrix

   for (int k=1; k<M; k++)

       // --- The matrix of Romberg approximations is triangular!
       for (int kk=0; kk<(M-k+1); kk++)
          // --- See the 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 the vector Rnum for numerical approximations
   thrust::host_vector<float> Rnum(M);
   thrust::copy(R.begin(), R.begin() + M, Rnum.begin());

   for (int i=0; i<M; i++) printf("%i %fn",i,Rnum[i]);

   return 0;

Leave a Reply

Your email address will not be published. Required fields are marked *