ArrayFire vs Thrust

test.cu

#include <stdio.h>
#include <arrayfire.h>
#include <thrust/device_vector.h>
#include <thrust/sort.h>
#include <thrust/binary_search.h>
#include <thrust/adjacent_difference.h>
using namespace af;
#define ITER 100

int main(int argc, char **argv)
{
    try {
        int num_particles = 2000;
        int num_cells = 1500;
        int *dev_particle_cell_indices = array::alloc<int>(num_particles);

        timer::tic();
        for (int i = 0; i < ITER; i++)
        {
            thrust::device_ptr<int> rawin = thrust::device_pointer_cast(dev_particle_cell_indices);

            // Sort a scratch copy of the cell indices by value
            thrust::device_vector<int> cidx(num_particles);
            thrust::copy(rawin, rawin+num_particles, cidx.begin());
            thrust::sort(cidx.begin(), cidx.end());

            // Use binary search to extract all the cell counts/offsets
            thrust::counting_iterator<int> cellnumber(0);
            thrust::device_vector<int> offsets(num_cells), counts(num_cells);

            // Offsets come from lower_bound of the ordered cell numbers
            thrust::lower_bound(cidx.begin(), cidx.end(), cellnumber, cellnumber+num_cells, offsets.begin());

            // Counts come from the adjacent_difference of the upper_bound of the ordered cell numbers
            thrust::upper_bound(cidx.begin(), cidx.end(), cellnumber, cellnumber+num_cells, counts.begin());
            thrust::adjacent_difference(counts.begin(), counts.end(), counts.begin());

        }
        printf("Thrust time taken: %lf\n",timer::toc()/ITER);

        timer::tic();
        for (int i = 0; i < ITER; i++)
        {
            array cell_indices(num_particles, 1, dev_particle_cell_indices, afDevicePointer);
            array particle_counts = zeros(num_cells);

            gfor(array i, num_cells) // Parallel for loop
                particle_counts(i) = sum(cell_indices == i);

            array particle_offsets = accum(particle_counts); // Inclusive sum
        }

        printf("ArrayFire time taken: %lf\n", timer::toc()/ITER);

    } catch (af::exception e) {
        std::cout << e.what();
    }

    return 0;
}

For those who are here and wondering about ArrayFire: It is free to use if you are connected to the internet.