アトミック型のヒストグラム集計アルゴリズムは、並行・並列処理に慣れている読者には馴染みの手法となります。

各ワークグループで256要素のビン（カテゴリ）を保持するローカルメモリの空間を割り当て、そこに複数のスレッドから同時アクセスが来ても整合性がとれるように、アトミック関数を使って処理をします。

サンプルコードではこれを以下のように記述しています。

for(int i = 0; i < BIN_SIZE; ++i)
{
   uchar value = data[group_id * local_size * BIN_SIZE + i * local_size + local_id];
   atomic_inc(&shared_local_memory[value]);
}

複数の並列するスレッド（ワークアイテム）がshared_local_memory変数をアクセスしますが、atomic_inc関数によってアトミックな加算処理が行われ正しい集計がとれます。

HistogramAtomicTest.py. 

import pyopencl as cl
import numpy as np
from numpy.random import *

# Set the seed to 100
np.random.seed(100)

INT_BYTES = 4
BIN_SIZE = 256
LOCAL_SIZE = 128
FACTOR = 64
DATA_SIZE = BIN_SIZE * LOCAL_SIZE * FACTOR
data = randint(0, 4, DATA_SIZE).astype(np.uint32)
out = np.zeros(BIN_SIZE*FACTOR).astype(np.uint32)
out.reshape((FACTOR, BIN_SIZE))

devices = [cl.get_platforms()[0].get_devices(cl.device_type.GPU)[0]]

ctx = cl.Context(devices)
queue = cl.CommandQueue(ctx)
mf = cl.mem_flags
data_mem = cl.Buffer(ctx, mf.USE_HOST_PTR, hostbuf=data)
out_mem = cl.Buffer(ctx, mf.WRITE_ONLY, size=out.nbytes)

program = cl.Program(ctx, """
    __kernel void histogram(
            __global const uint* data,
            __local uint* shared_local_memory,
            __global uint* buckets)
    {

         size_t local_id   = get_local_id(0);
         size_t group_id   = get_group_id(0);
         size_t local_size = get_local_size(0);

         shared_local_memory[local_id] = 0;
         uint4* input = (uint4*) shared_local_memory;
         for(int i = 0; i < 64; i++)
             input[i] = 0;
         barrier(CLK_LOCAL_MEM_FENCE);

         for(int i = 0; i < BIN_SIZE; ++i)
         {
           uchar value = data[group_id * local_size * BIN_SIZE + i * local_size + local_id];
           atomic_inc(&shared_local_memory[value]);
         }
         barrier(CLK_LOCAL_MEM_FENCE);

         // 0...127
         if(local_id == 0) {
             for(int i = 0; i < BIN_SIZE; i++) {
               buckets[group_id * BIN_SIZE  + i] = shared_local_memory[i]; // 0..63 * 256 + 0...255
             }
         }
    }
    """).build(options=["-DBIN_SIZE=256"])

kernel = cl.Kernel(program, name="histogram")

kernel.set_arg(0, data_mem)
# size to be specified in bytes (ie., 32 bits -> 4 bytes)
kernel.set_arg(1, cl.LocalMemory(INT_BYTES*BIN_SIZE))
kernel.set_arg(2, out_mem)

cl.enqueue_nd_range_kernel(
    queue=queue,
    kernel=kernel,
    global_work_size=(FACTOR* LOCAL_SIZE, 1, 1),
    local_work_size=(LOCAL_SIZE, 1, 1))

out = cl.enqueue_map_buffer(
    queue=queue,
    buf=out_mem,
    flags=mf.WRITE_ONLY,
    offset=0,
    shape=(FACTOR, BIN_SIZE,),
    dtype=np.uint32
)

#print(out[0][0][0:256])
print(sum(sum(out[0])))

上記のプログラムの出力は以下のようになります。

出力. 

2097152

ヒストグラムの総和が一致したことが確認できます。