A faster way to compute an array in thrust-cuda

Question

I have four arrays:

thrust::device_vector<int> vertices;
thrust::device_vector<int> adjacency:
thrust::device_vector<int> degree_summation;
thrust::device_vector<int> degree;

The indexes of vertices array represent the vertices of a graph while its content points to the first neighbour of each vertex in the adjacency array. The array degree stores the degree of each vertex in vertices.

I want the degree_summation array keeps the degree of each vertex plus the summation of the degree of each vertex neighbourhood.

Ex: Given the graph:
 0
/ \
1--2

vertices = {0,2,4}
adjacency = {1,2,0,2,0,1}
degree = {2,2,2}

What i want:

degree_summation = {6,6,6}

Currently, i am using a for loop to compute this value, but i think perhaps there is a faster way to compute these values using the primitives given by thrust. The kernel i use to compute the degree_sumation array:

__global__ void computeDegreeSummation(int* vertices,int* adjacency,unsigned int* degree_summation, unsigned int * degree,unsigned int num_vertices, unsigned int num_edges){
unsigned long int tid = (gridDim.y*blockIdx.y+blockIdx.x)*blockDim.x+threadIdx.x;
if(tid < num_vertices){
    int pos_first_neighbor = vertices[tid];
    int pos_last_neighbor;
    if(tid != num_vertices - 1){
        int it = tid;
        pos_last_neighbor = vertices[it+1];
        while(pos_last_neighbor == 0){
            it++;
            pos_last_neighbor = vertices[it+1];
        }
        pos_last_neighbor--;
    }//if
    else{
        pos_last_neighbor = num_edges - 1;
    }//else
    for(int nid = pos_first_neighbor; nid <= pos_last_neighbor; nid++){
        if(adjacency[nid]!=tid){
            degree_summation[tid]+=degrees[adjacency[nid]];
        }//if
    }//for
}//if
}//kernel

There is a if inside the for since i had initialized the degree_summation with the degree of each vertex before compute the kernel.

Answer 1

I don't know if this will be faster - since I don't know how fast your reference implementation runs. But this is one possible approach:

1. create a vector which defines the neighborhoods in adjacency for each vertex.

   vertices =     {0,  2,  4}
   adjacency =    {1,2,0,2,0,1}
   neighborhood = {1,0,1,0,1,0}  (after scatter)
   neighborhood = {1,1,2,2,3,3}  (after inclusive scan)
   

   To accomplish this, we will thrust::scatter a set of 1's (thrust::constant_iterator) into a neighborhood array according to the indices provided by vertices, then do a thrust::inclusive_scan on neighborhood.

2. We can then use thrust::reduce_by_key using the neighborhood array as the keys, and using a thrust::permutation_iterator to select the value out of degree corresponding to each vertex listed in each neighborhood of adjacency.

3. We then just need to add degree to the intermediate degree_summation result produced in step 2. We can use thrust::transform for this.

Here is a fully worked example:

$ cat t612.cu
#include <thrust/device_vector.h>
#include <thrust/iterator/constant_iterator.h>
#include <thrust/scatter.h>
#include <thrust/scan.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/reduce.h>
#include <thrust/iterator/discard_iterator.h>
#include <thrust/transform.h>
#include <thrust/functional.h>
#include <thrust/copy.h>
#include <iostream>


int main(){

  // input data setup

  int h_vertices[] = {0,2,4};
  int h_adjacency[] = {1,2,0,2,0,1};
  int h_degree[] = {2,2,2};
  int vertices_size = sizeof(h_vertices)/sizeof(int);
  int adjacency_size = sizeof(h_adjacency)/sizeof(int);
  int degree_size = sizeof(h_degree)/sizeof(int);
  thrust::device_vector<int> vertices(h_vertices, h_vertices+vertices_size);
  thrust::device_vector<int> adjacency(h_adjacency, h_adjacency+adjacency_size);
  thrust::device_vector<int> degree(h_degree, h_degree+degree_size);
  thrust::device_vector<int> degree_summation(degree_size);

  // create neighborhood array

  thrust::device_vector<int> neighborhood(adjacency_size);
  thrust::constant_iterator<int> first(1);
  thrust::scatter(first, first+adjacency_size, vertices.begin(), neighborhood.begin());
  thrust::inclusive_scan(neighborhood.begin(), neighborhood.end(), neighborhood.begin());

  // sum degree over the neighborhoods

  thrust::reduce_by_key(neighborhood.begin(), neighborhood.end(), thrust::make_permutation_iterator(degree.begin(), adjacency.begin()), thrust::make_discard_iterator(), degree_summation.begin());

  // add the vertex degrees

  thrust::transform(degree.begin(), degree.end(), degree_summation.begin(), degree_summation.begin(), thrust::plus<int>());

  // display results

  thrust::copy(degree_summation.begin(), degree_summation.end(), std::ostream_iterator<int>(std::cout, ","));
  std::cout << std::endl;
  return 0;
}
$ nvcc -arch=sm_35 -o t612 t612.cu
$ ./t612
6,6,6,
$

(deleting the rest of my previous error. the alternative approach I had suggested was not valid.)