用户
 找回密码
 立即注册
91金属狂人 该用户已被删除
发表于 2013-10-21 11:08:09
69797
本来想学习一下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]++;
    }
}



使用道具 举报 回复
发表于 2013-10-21 11:44:16
本帖最后由 横扫千军 于 2013-10-21 11:46 编辑

楼主您好,

的确,您给出的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的保证,这个无视即可。

感谢您的来访
使用道具 举报 回复 支持 反对
发表于 2013-10-21 13:40:28
横扫千军 发表于 2013-10-21 11:44
楼主您好,

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

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

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

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

感谢来访。
使用道具 举报 回复 支持 反对
发表于 2013-10-21 14:06:21
横扫千军 发表于 2013-10-21 14:03
恭喜楼主学会了大部分内容。

关于为何一个1B的类型,无视最低2bit后,就是将0-3为一类,4-7为一类.....

恩,版主这样一说,我就明白了,谢谢了。我往下看CUDA版本了,若有不明白,再来请教。
使用道具 举报 回复 支持 反对
发表于 2013-10-22 10:10:43
91金属狂人 发表于 2013-10-21 14:06
恩,版主这样一说,我就明白了,谢谢了。我往下看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;
    }
}
使用道具 举报 回复 支持 反对
发表于 2013-10-22 10:54:33
横扫千军 发表于 2013-10-21 14:03
恭喜楼主学会了大部分内容。

关于为何一个1B的类型,无视最低2bit后,就是将0-3为一类,4-7为一类.....

    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,最后0*s_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类。                                                
使用道具 举报 回复 支持 反对
发表于 2013-10-22 11:32:39
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;
使用道具 举报 回复 支持 反对
发新帖
您需要登录后才可以回帖 登录 | 立即注册