Using Thrust library of CUDA for large values

Question

Hi I wanted to implement a loop which is extremely large in thrust but i find it much slower than normal C++ code. Can you please tell me where am i going wrong. fi and fj are host vectors

xsize usually is around a 7-8 digit number

thrust::host_vector <double> df((2*floor(r)*(floor(r)+1)+1)*n*n); 
thrust::device_vector<double> gpu_df((2*floor(r)*(floor(r)+1)+1)*n*n); 
     for(i=0;i<xsize;i++)
     {
        gpu_df[i]=(fi[i]-fj[i]);

         if(gpu_df[i]<0)
            gpu_df[i]=0;
        else 
       gpu_df[i]=gpu_df[i]*(fi[i]-fj[i]);
        if(gpu_df[i]>255)
            gpu_df[i]=255;
        //      cout<<fi[i]<<"\n";
     }
df=gpu_df;

I feel the code is not being parallelized. Could you please help me out.

Answer 1

To run programs on the GPU with Thrust you need to write them in terms of Thrust algorithms like reduce, transform, sort, etc. In this case we can write the computation in terms of transform, since the loop was just computing a function F(fi[i], fj[i]) and storing the result in df[i]. Note that we must first move the input arrays to the device before calling transform because Thrust requires the input and output arrays to live in the same place.

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/functional.h>
#include <cstdio>

struct my_functor
  : public thrust::binary_function<float,float,float>
{
  __host__ __device__
  float operator()(float fi, float fj)
      {
    float d =  fi - fj;

    if (d < 0)
      d = 0;
    else
      d = d * d;

    if (d > 255)
      d = 255;

    return d;
  }
};

int main(void)
{
  size_t N = 5;

  // allocate storage on host
  thrust::host_vector<float>   cpu_fi(N);
  thrust::host_vector<float>   cpu_fj(N);
  thrust::host_vector<float>   cpu_df(N);

  // initialze fi and fj arrays
  cpu_fi[0] = 2.0;  cpu_fj[0] =  0.0;
  cpu_fi[1] = 0.0;  cpu_fj[1] =  2.0;
  cpu_fi[2] = 3.0;  cpu_fj[2] =  1.0;
  cpu_fi[3] = 4.0;  cpu_fj[3] =  5.0;
  cpu_fi[4] = 8.0;  cpu_fj[4] = -8.0;

  // copy fi and fj to device
  thrust::device_vector<float> gpu_fi = cpu_fi;
  thrust::device_vector<float> gpu_fj = cpu_fj;

  // allocate storage for df
  thrust::device_vector<float> gpu_df(N);

  // perform transformation
  thrust::transform(gpu_fi.begin(), gpu_fi.end(),  // first input range
                    gpu_fj.begin(),                // second input range
                    gpu_df.begin(),                // output range
                    my_functor());                 // functor to apply

  // copy results back to host
  thrust::copy(gpu_df.begin(), gpu_df.end(), cpu_df.begin());

  // print results on host
  for (size_t i = 0; i < N; i++)
    printf("f(%2.0lf,%2.0lf) = %3.0lf\n", cpu_fi[i], cpu_fj[i], cpu_df[i]);

  return 0;
}

For reference, here's the output of the program:

f( 2, 0) =   4
f( 0, 2) =   0
f( 3, 1) =   4
f( 4, 5) =   0
f( 8,-8) = 255