CUDA histogram reduce_by_key failing

Question

I have the following CUDA Thrust code that uses reduce_by_key to histogram the values [0, 1024) into 256 buckets. I expect each bucket to have a count = 4, yet I see bucket 0 has 256, bucket 255 has 3, and the remainder have 4.

#include <stdio.h>
#include <stdlib.h>

#include <cuda.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>

#include <thrust/device_vector.h>
#include <thrust/extrema.h>
#include <thrust/pair.h>

#define SIZE 1024

struct binFunc {
    const float minVal;
    const float valRange;
    const int numBins;
    binFunc(float _minVal, float _valRange, int _numBins) :
        minVal(_minVal), valRange(_valRange), numBins(_numBins) {}

    __host__ __device__
    int operator()(float v) const {
        int b = int((v - minVal) / valRange * float(numBins));
        return b;
    }
};

int main() {
    thrust::device_vector<float> d_vec(SIZE);
    for (int i = 0; i < SIZE; ++i)
        d_vec[i] = float(i);

    thrust::device_vector<float>::iterator min;
    thrust::device_vector<float>::iterator max;
    thrust::pair<thrust::device_vector<float>::iterator,
            thrust::device_vector<float>::iterator> minmax =
            thrust::minmax_element(d_vec.begin(), d_vec.end());
    min = minmax.first;
    max = minmax.second;
    float minVal = *min;
    float maxVal = *max;

    std::cout << "The minimum value is " << minVal
            << " and the maximum value is " << maxVal << "." << std::endl;

    float valRange = maxVal - minVal;

    std::cout << "The range is " << valRange << "." << std::endl;

    int numBins = 256;

    thrust::device_vector<int> d_binResults(SIZE);
    thrust::transform(d_vec.begin(), d_vec.end(), d_binResults.begin(),
            binFunc(minVal, valRange, numBins));

    thrust::device_vector<int>::iterator d_binResults_iter =
            d_binResults.begin();
    for (int i = 0; i < 10; ++i) {
        int b = *d_binResults_iter;
        printf("d_binResults[%d]=%d\n", i, b);
        d_binResults_iter++;
    }

    std::cout << "The numBins is " << numBins << "." << std::endl;

    thrust::device_vector<int> d_binsKeys(numBins);
    thrust::device_vector<int> d_binsValues(numBins);

    thrust::pair<thrust::device_vector<int>::iterator,
            thrust::device_vector<int>::iterator> keys_and_values =
            thrust::reduce_by_key(d_binResults.begin(), d_binResults.end(),
                    thrust::constant_iterator<int>(1), d_binsKeys.begin(),
                    d_binsValues.begin());

    thrust::device_vector<int>::iterator d_binsKeys_begin_iter =
            d_binsKeys.begin();
    thrust::device_vector<int>::iterator d_binsValues_begin_iter =
            d_binsValues.begin();
    for (int i = 0; i < numBins; ++i) {
        int key = *d_binsKeys_begin_iter;
        int val = *d_binsValues_begin_iter;
        printf("d_binsValues[%d]=(%d,%d)\n", i, key, val);
        d_binsKeys_begin_iter++;
        d_binsValues_begin_iter++;
    }
    return 0;
}

The salient part of the output is:

d_binsValues[0]=(0,256)
d_binsValues[1]=(1,4)
d_binsValues[2]=(2,4)
...
d_binsValues[254]=(254,4)
d_binsValues[255]=(255,3)

So, bucket 0 has 256 elements, and bucket 255 has 3 elements? What's going on here?

Answer 1

If you print out all the d_binResults[] values instead of the first 10, you will discover that the last element (d_binResults[1023]) has a value of 256! But that is an invalid bin index. For numBins = 256, the valid indices are 0..255.

It is occurring due to the calculation arithmetic in your functor:

    int b = int((v - minVal) / valRange * float(numBins));

Plugging in the relevant values for the last element, we have:

(1023 - 0)/1023*256 = 256

But 256 is an invalid bin index. It turns out that this breaks the reduce_by_key operation, causing both the last bin to have 3 elements and the first bin to be "corrupted".

If you fix this you will fix both issues you describe (first bin has 256 elements, last bin has 3.)

As a simple proof, add this line of code:

d_binResults[1023] = 255;

immediately after your thrust::transform operation. The results are then correct. How you choose to correct your bin calculation arithmetic is up to you. (possibly "fixable" by adding 1 to valRange but that may imply something about your expected histogram values).