想学习一下SDK里的histogram例子

本来想学习一下CUDA里的histogram,最简单的CPU写法就是,for(int i = 0; i < BIN_COUNT; i++) result[i] = 0;
for(int i = 0; i < dataN; i++)
result[data[i]]++;
但histogram的CPU版本看不明白。其中这一句 assert(sizeof(uint) == 4 && (byteCount % 4) == 0);不明白,还有为什么要把输入的字节类型数据转换为无符号整形,uint data = ((uint *)h_Data)[i];
这4行是进行字节操作吗?
h_Histogram[(data >> 2) & 0x3FU]++; h_Histogram[(data >> 10) & 0x3FU]++;
h_Histogram[(data >> 18) & 0x3FU]++;
h_Histogram[(data >> 26) & 0x3FU]++;

代码如下,请版主指教一二,
详细代码可以在 C:\ProgramData\NVIDIA Corporation\CUDA Samples\v5.5\3_Imaging\histogram 找到

extern “C” void histogram64CPU(
uint *h_Histogram,
void *h_Data,
uint byteCount
)
{
for (uint i = 0; i < HISTOGRAM64_BIN_COUNT; i++)
h_Histogram[i] = 0;

assert(sizeof(uint) == 4 && (byteCount % 4) == 0);

for (uint i = 0; i < (byteCount / 4); i++)
{
uint data = ((uint *)h_Data)[i];
h_Histogram[(data >> 2) & 0x3FU]++;
h_Histogram[(data >> 10) & 0x3FU]++;
h_Histogram[(data >> 18) & 0x3FU]++;
h_Histogram[(data >> 26) & 0x3FU]++;
}
}

楼主您好,

的确,您给出的result[data[ i ]]++; 也是常用的分类形式一种。当data[ i ]的范围是[0,255]的时候,可以给出256种分类。

但是楼主您要知道,并不是所有的算法都强制要求必须分成256种的。分类的总数量被实际需要所决定,而不是固定的256种。
而数据类型,则要看原始的数据类型是什么,以及,访问的策略如何,也同样不总是固定成char的。

回到楼主您的问题上,则我们可以从代码中看出:
(1)总是使用uint读取,可能是为了将4个1B的类型打包,减少访存的次数。
(2)使用unsigned则为了防止编译器进行算术右移,防止进行错误的符号扩展。
(3)为何是:
h_Histogram[(data >> 2) & 0x3FU]++;
h_Histogram[(data >> 10) & 0x3FU]++;
h_Histogram[(data >> 18) & 0x3FU]++;
h_Histogram[(data >> 26) & 0x3FU]++;
可能该算法就一共需要64种类别。

对于里面的每个4B的值中的4个1B,
即:数据的位[2, 8), [10, 16), [18, 24), [26, 32)。
进行4次统计,每次只需要高6位,忽略低2位。

诚然,如果我们统计了[0,8), [9,16), [17, 24), [25,32)
(也就是你的直接用的byte的数组,然后hist[data[ i ]] ++;)
的确也可以。

但上文说了,一是人家只需要64种,你不能说你用过256种的,就强制人家不能用64种的。
而使用64种的,舍弃低有效的2位值,也是人之常情。(8位有效数据,如果你需要6位,你总不能舍弃最高有效的2位吧。那样就成错误了)。

二是前文说了,一次打包读取4个,方便简单高效。

然后回到楼主的最后一个小问题:
(4)为何加入assert里的2个sizeof(uint) == 4, 和byteCount % 4 == 0这2个条件,
则是后者可能是算法约定好的数据总量是4的倍数,这里强制验证下。
而前者,sizeof(uint) == 4,则是不好的代码风格,当使用者不确定当前编译器的uint是否是固定的32位值的时候,他应该使用uint32_t来替代(需要#include<stdint.h>,而不是恳请编译器的怜悯和assert的保证,这个无视即可。

感谢您的来访
[/i]

恩,现在就是要把数组h_Data(类型为uchar,数据在0-255之间)分为64类,即0到3为一类,4至7为一类… …252至255为一类,也即是相邻4个数为一类。把输入的字节类型数据转换为无符号整形,是为了将4个1B的类型打包,减少访存的次数,然后再用后4行来进行每一个B分类。前面都明白了。数据右移2位再与64,就是用来分成64类。我手动移了十几次,发现做出来的结果跟计算机做出来的一样的,我明白了最后4句的用途。但是还没能理解。谢谢版主指引了。我会常来的,呵呵。

恭喜楼主学会了大部分内容。

关于为何一个1B的类型,无视最低2bit后,就是将0-3为一类,4-7为一类…
这其实很简单,最后的2个bit能表示4个数。自然无视后就每4个数的类别是一样的。
你别看它是十进制的数,当成2进制后,就能立刻明白了。

感谢来访。

恩,版主这样一说,我就明白了,谢谢了。我往下看CUDA版本了,若有不明白,再来请教。

版主,你好,昨天看了半天,只能理解CUDA版本的部分代码,请版主讲解一下思路。
版主应该对SDK里的大多数例子都挺熟悉吧?这是SDK里的例子
详细代码可以在 C:\ProgramData\NVIDIA Corporation\CUDA Samples\v5.5\3_Imaging\histogram
//Data type used for input data fetches
typedef uint4 data_t;

//May change on future hardware, so better parametrize the code
define SHARED_MEMORY_BANKS 16

////////////////////////////////////////////////////////////////////////////////
// Main computation pass: compute gridDim.x partial histograms
////////////////////////////////////////////////////////////////////////////////
//Count a byte into shared-memory storage
inline device void addByte(uchar *s_ThreadBase, uint data)
{
s_ThreadBase[UMUL(data, HISTOGRAM64_THREADBLOCK_SIZE)]++;
}

//Count four bytes of a word
inline device void addWord(uchar *s_ThreadBase, uint data)
{
//Only higher 6 bits of each byte matter, as this is a 64-bin histogram
addByte(s_ThreadBase, (data >> 2) & 0x3FU); 这里跟CPU版本的相似,找到第几类再乘线程块
addByte(s_ThreadBase, (data >> 10) & 0x3FU); 宽度64的那个位置加1,为什么是64,难道是把
addByte(s_ThreadBase, (data >> 18) & 0x3FU); 类数存在一列里。
addByte(s_ThreadBase, (data >> 26) & 0x3FU);
}

global void histogram64Kernel(uint *d_PartialHistograms, data_t *d_Data, uint dataCount)
{
//Encode thread index in order to avoid bank conflicts in s_Hist access:
//each group of SHARED_MEMORY_BANKS threads accesses consecutive shared memory banks
//and the same bytes [0…3] within the banks
//Because of this permutation block size should be a multiple of 4 * SHARED_MEMORY_BANKS
const uint threadPos =
((threadIdx.x & ~(SHARED_MEMORY_BANKS * 4 - 1)) << 0) | 每个线程对应不同的threadPos
((threadIdx.x & (SHARED_MEMORY_BANKS - 1)) << 2) | 线程0对应threadPos=0,线程1
((threadIdx.x & (SHARED_MEMORY_BANKS * 3)) >> 4); 对应threadPos=4,线程2对应
threadPos=8,应该是避免bank冲突
//Per-thread histogram storage
shared uchar s_Hist[HISTOGRAM64_THREADBLOCK_SIZE * HISTOGRAM64_BIN_COUNT];
uchar *s_ThreadBase = s_Hist + threadPos;
声明为字节型共享内存,并且对应的地址有对应的threadPos。
//Initialize shared memory (writing 32-bit words)
#pragma unroll

for (uint i = 0; i < (HISTOGRAM64_BIN_COUNT / 4); i++)
{
((uint *)s_Hist)[threadIdx.x + i * HISTOGRAM64_THREADBLOCK_SIZE] = 0;
} 将声明的字节型共享内存转为uint,用4字来初始化,当i=0时,threadIdx.x 0-63,初始化了64个4字了。

//Read data from global memory and submit to the shared-memory histogram
//Since histogram counters are byte-sized, every single thread can’t do more than 255 submission
__syncthreads();

for (uint pos = UMAD(blockIdx.x, blockDim.x, threadIdx.x); pos < dataCount; pos += UMUL(blockDim.x, gridDim.x))
{ 这里应该是读数据进共享内存了,但不怎么明白。
data_t data = d_Data[pos];
addWord(s_ThreadBase, data.x);
addWord(s_ThreadBase, data.y);
addWord(s_ThreadBase, data.z);
addWord(s_ThreadBase, data.w);
}
这下面应该是将共享内存里的数据进行统计吧?
//Accumulate per-thread histograms into per-block and write to global memory
__syncthreads();

if (threadIdx.x < HISTOGRAM64_BIN_COUNT)
{
uchar *s_HistBase = s_Hist + UMUL(threadIdx.x, HISTOGRAM64_THREADBLOCK_SIZE);

uint sum = 0;
uint pos = 4 * (threadIdx.x & (SHARED_MEMORY_BANKS - 1));

#pragma unroll

for (uint i = 0; i < (HISTOGRAM64_THREADBLOCK_SIZE / 4); i++)
{
sum +=
s_HistBase[pos + 0] +
s_HistBase[pos + 1] +
s_HistBase[pos + 2] +
s_HistBase[pos + 3];
pos = (pos + 4) & (HISTOGRAM64_THREADBLOCK_SIZE - 1);
}

d_PartialHistograms[blockIdx.x * HISTOGRAM64_BIN_COUNT + threadIdx.x] = sum;
}
}

for (uint pos = UMAD(blockIdx.x, blockDim.x, threadIdx.x); pos < dataCount; pos += UMUL(blockDim.x, gridDim.x))
{
data_t data = d_Data[pos]; 不知道这样理解对否,当pos=0时,读得16字,当执行第一个
addWord(s_ThreadBase, data.x); addWord时,假设data.x由(data >> 2) & 0x3FU 算出是第0
addWord(s_ThreadBase, data.y); 类,然后0乘64,当pos=0时,threadPos=0,也即是
addWord(s_ThreadBase, data.z); s_ThreadBase=s_Hist+0,最后0s_ThreadBase位置加1
addWord(s_ThreadBase, data.w); 当pos=0时,threadPos=4,也即是s_ThreadBase=s_Hist+4
} 假设data.x由(data >> 2) & 0x3FU 算出是第1哦的,然后1乘64 ,对应1
s_ThreadBase位置
加1,这样保证了一个warp内线程不会访问同一个bank,但又能统计64类。共享内存为64行*64列,其中行应该代表类别,第0行代表第1类,第1行,代表第1类。

if (threadIdx.x < HISTOGRAM64_BIN_COUNT)
{
uchar *s_HistBase = s_Hist + UMUL(threadIdx.x, HISTOGRAM64_THREADBLOCK_SIZE);

uint sum = 0;
uint pos = 4 * (threadIdx.x & (SHARED_MEMORY_BANKS - 1));

#pragma unroll
这里对第行进行求和,线程0对第一行进行求各,线程1对第二行进行求和,一次求和求4个字,故for循环16次,在一个warp内,通过uint pos = 4 * (threadIdx.x & (SHARED_MEMORY_BANKS - 1)); 使得不会bank冲突,求得的第0行,即第0类,第1行,即第1类,分别把sum写回 d_PartialHistograms。 总算明白了这个核函数,呵呵,版主还没回答,我就自己想明白了。有不明白的地方,我会再来的,谢谢版主。

for (uint i = 0; i < (HISTOGRAM64_THREADBLOCK_SIZE / 4); i++)
{
sum +=
s_HistBase[pos + 0] +
s_HistBase[pos + 1] +
s_HistBase[pos + 2] +
s_HistBase[pos + 3];
pos = (pos + 4) & (HISTOGRAM64_THREADBLOCK_SIZE - 1);
}

d_PartialHistograms[blockIdx.x * HISTOGRAM64_BIN_COUNT + threadIdx.x] = sum;