Nézzük először meg a program keretét. A jelen program egy véletlenszerű adatsort generál, aminek utána kiszámítja a hisztogramját.

Az adatok előkészítése a bemeneti kép és a hisztogram tömb előkészítéséből áll.

    unsigned char *buffer =
                     (unsigned char*)big_random_block( SIZE );

    unsigned char *dev_buffer;
    unsigned int *dev_histo;
    cudaMalloc( (void**)&dev_buffer, SIZE );
    cudaMemcpy( dev_buffer, buffer, SIZE, cudaMemcpyHostToDevice );

    cudaMalloc( (void**)&dev_histo, 256 * sizeof( int ) );
    cudaMemset( dev_histo, 0, 256 * sizeof( int ) );

Az adat előkészítés után már jöhet is a kernel indítása. A kernelt megint csak belső szorzat számításához hasonlóan indítjuk őgy, hogy egy szál több adatot is feldolgozhat.

    cudaDeviceProp  prop;
    cudaGetDeviceProperties( &prop, 0 );
    int blocks = prop.multiProcessorCount;
    histo_kernel<<<blocks*2,256>>>( dev_buffer, SIZE, dev_histo );

A program most annyiban van "megbolondítva", hogy először lekérdezi a multiprocesszorok számát, és kétszer ennyi blokkot indít, blokkonként 256 szállal. Ezzel a futás közben dinamikusan tudjuk a GPU-hoz igazítani az indított szálstruktúrát.

A program keret többi része csak utómunkákat ír le térjünk inkább rá, hogy mit is csinál a kernel.

__global__ void histo_kernel( unsigned char *buffer,
                              long size,
                              unsigned int *histo ) {
    // calculate the starting index and the offset to the next
    // block that each thread will be processing
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    int stride = blockDim.x * gridDim.x;
    while (i < size) {
        atomicAdd( &histo[buffer[i]], 1 );
        i += stride;
    }
}

A kernel szintén nem túl bonyolult. A belsőszorzathoz hasonlóan a szálak együttesen végigolvassák az adatsort, és minden adatelemnél (vagy pixelnél) nivelik a hisztogram megfelelő számlálóját egyet. A növelés az atomicAdd(...)-al történik, így a végső szám garantáltan helyes lesz.

Ezen a ponton meg kell viszont jegyeznünk, hogy az atomikus műveletek lassítják a programunkat. A GPU nagy teljesítmény a számítások párhuzamos végzésében áll. Egy memória elérés kb. 600 órajelet vesz igénybe (egy írás-olvasásá pár pedig 1200-at), így ha a szálak atomikus műveletek miatt egymásra várnak, az nagyban rontja a hatákonyságot. Ezért atomikus műveleteket csak akkor használjunk, ha tényleg elkerülhetetlen.

/* This software contains source code provided by NVIDIA Corporation.
 * The program is based on the example code hist_gpu_gmem_atomics.cu of
 * Chapter 9 of the "CUDA by Example" book.
 * 
 * A program az NVIDIA Corporation által készített programkódot tartalmaz.
 * A program a "Cuda by Example" tankönyv 9. hist_gpu_gmem_atomics.cu 
 * példaprogramjára épül.
 */

#include 
#include 

#include 
#include 
#include 

using namespace std;

#define SIZE    (100*1024*1024)

void* big_random_block( int size ) {
    unsigned char *data = (unsigned char*)malloc( size );
    for (int i=0; i<size; i++)
        data[i] = rand()%256;

    return data;
}

__global__ void histo_kernel( unsigned char *buffer,
                              long size,
                              unsigned int *histo ) {
    // calculate the starting index and the offset to the next
    // block that each thread will be processing
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    int stride = blockDim.x * gridDim.x;
    while (i < size) {
        atomicAdd( &histo[buffer[i]], 1 );
        i += stride;
    }
}

int main( void ) {
    unsigned char *buffer =
                     (unsigned char*)big_random_block( SIZE );

    // capture the start time
    // starting the timer here so that we include the cost of
    // all of the operations on the GPU.
    cudaEvent_t     start, stop;
    cudaEventCreate( &start );
    cudaEventCreate( &stop );
    cudaEventRecord( start, 0 );

    // allocate memory on the GPU for the file's data
    unsigned char *dev_buffer;
    unsigned int *dev_histo;
    cudaMalloc( (void**)&dev_buffer, SIZE );
    cudaMemcpy( dev_buffer, buffer, SIZE, cudaMemcpyHostToDevice );

    cudaMalloc( (void**)&dev_histo, 256 * sizeof( int ) );
    cudaMemset( dev_histo, 0, 256 * sizeof( int ) );

    // kernel launch - 2x the number of mps gave best timing
    cudaDeviceProp  prop;
    cudaGetDeviceProperties( &prop, 0 );
    int blocks = prop.multiProcessorCount;
    histo_kernel<<<blocks*2,256>>>( dev_buffer, SIZE, dev_histo );
    
    unsigned int    histo[256];
    cudaMemcpy( histo, dev_histo, 256 * sizeof( int ), cudaMemcpyDeviceToHost );

    // get stop time, and display the timing results
    cudaEventRecord( stop, 0 );
    cudaEventSynchronize( stop );
    float   elapsedTime;
    cudaEventElapsedTime( &elapsedTime, start, stop );

    long histoCount = 0;
    for (int i=0; i<256; i++) {
        histoCount += histo[i];
        cout << i << ": " << histo[i] << endl;
    }
    
    cout << endl;
    cout << "Time to generate:  " << elapsedTime << " msec." << endl;
    cout << "Histogram Sum: " << histoCount << endl;

    // verify that we have the same counts via CPU
    for (int i=0; i<SIZE; i++)
        histo[buffer[i]]--;
    for (int i=0; i<256; i++) {
        if (histo[i] != 0)
            cout << "Failure at " << i << "!  Off by " << histo[i] << endl;
    }

    cudaEventDestroy( start );
    cudaEventDestroy( stop );
    cudaFree( dev_histo );
    cudaFree( dev_buffer );
    free( buffer );
    return 0;
}