感興趣的朋友可以從http://code.google.com/p/imagefilter-opencl/downloads/detail?name=amdunicourseCode2.zip&can=2&q=#makechanges下載完整代碼。

注意代碼運行后，會在程序目錄生成lenna_rotate.jpg,這時gpu執(zhí)行的結(jié)果，另外還有一個cpu_lenna_rotate.jpg這是CPU執(zhí)行的結(jié)果。

3、一個矩陣乘法的例子

??? 在amd的slides中，本節(jié)還講了一個簡單的，沒有優(yōu)化的矩陣乘法，一共才2兩頁ppt，所以我也不在這兒詳細(xì)講述了，…,但簡單介紹還是需要的。

上面的代碼是矩陣乘法的例子，有三重循環(huán)，下面我們只給出kernel代碼，完整程序請從：http://code.google.com/p/imagefilter-opencl/downloads/detail?name=amdunicodeCode3.zip&can=2&q=#makechanges下載。

從今天開始學(xué)習(xí)OpenCL……

??????? 因為老狼的顯卡是AMD 5xx的redwood,所以下面先介紹OpenCL APP(Accelerated Parallel processing)的安裝。

下載地址：http://developer.amd.com/tools/hc/AMDAPPSDK/downloads/Pages/default.aspx

安裝注意事項：http://developer.amd.com/tools/hc/AMDAPPSDK/assets/AMD_APP_SDK_Installation_Notes.pdf

有用的資料：http://developer.amd.com/tools/hc/AMDAPPSDK/documentation/Pages/default.aspx

其中最有用的是下面幾個文檔：

???? AMD最新顯卡Tahiti的ISA介紹，對OpenCL編程優(yōu)化有用。

http://developer.amd.com/wordpress/media/2012/12/AMD_Southern_Islands_Instruction_Set_Architecture.pdf

?

????? 這本書是最好的OpenCL教程，好過市面上的任何一本OpenCL書，其中包括很多優(yōu)化Kernel代碼的內(nèi)容，我計劃以后就按照這本書的內(nèi)容來學(xué)習(xí)opencl。

http://developer.amd.com/download/AMD_Accelerated_Parallel_Processing_OpenCL_Programming_Guide.pdf

?

????? 再就是OpenCL 1.2的spec了，下載地址：OpenCL? 1.2 Specification (revision 15) ，相對于1.1來說，1.2中還是有一些變化的，比如我以前寫的程序中CreateImage2D函數(shù)發(fā)現(xiàn)在1.2中沒有了，spec其實就是一個函數(shù)手冊，偶爾用來查詢一下而已。

??? 另外一個比較好的初級教程，就是我翻譯的AMD OpenCL大學(xué)教程了，http://www.cnblogs.com/mikewolf2002/archive/2012/01/30/2332356.html

現(xiàn)在，我們開始寫一個簡單的OpenCL程序，計算兩個數(shù)組相加的和，放到另一個數(shù)組中去。程序用cpu和gpu分別計算，最后驗證它們是否相等。OpenCL程序的流程大致如下：

下面是source code中的主要代碼：

?

int main(int argc, char* argv[])
??? {
??? //在host內(nèi)存中創(chuàng)建三個緩沖區(qū)
??? float *buf1 = 0;
??? float *buf2 = 0;
??? float *buf = 0;

??? buf1 =(float *)malloc(BUFSIZE * sizeof(float));
??? buf2 =(float *)malloc(BUFSIZE * sizeof(float));
??? buf =(float *)malloc(BUFSIZE * sizeof(float));

??? //用一些隨機(jī)值初始化buf1和buf2的內(nèi)容
??? int i;
??? srand( (unsigned)time( NULL ) );
??? for(i = 0; i < BUFSIZE; i++)
??????? buf1[i] = rand()%65535;

??? srand( (unsigned)time( NULL ) +1000);
??? for(i = 0; i < BUFSIZE; i++)
??????? buf2[i] = rand()%65535;

??? //cpu計算buf1,buf2的和
??? for(i = 0; i < BUFSIZE; i++)
??????? buf[i] = buf1[i] + buf2[i];

??? cl_uint status;
??? cl_platform_id platform;

??? //創(chuàng)建平臺對象
??? status = clGetPlatformIDs( 1, &platform, NULL );

????? 注意：如果我們系統(tǒng)中安裝不止一個opencl平臺，比如我的os中，有intel和amd兩家opencl平臺，用上面這行代碼，有可能會出錯，因為它得到了intel的opencl平臺，而intel的平臺只支持cpu，而我們后面的操作都是基于gpu，這時我們可以用下面的代碼，得到AMD的opencl平臺。
?

??? //創(chuàng)建GPU設(shè)備
??? clGetDeviceIDs( platform, CL_DEVICE_TYPE_GPU,
??????? 1,
??????? &device,
??????? NULL);
??? //創(chuàng)建context
??? cl_context context = clCreateContext( NULL,
??????? 1,
??????? &device,
??????? NULL, NULL, NULL);
??? //創(chuàng)建命令隊列
??? cl_command_queue queue = clCreateCommandQueue( context,
??????? device,
??????? CL_QUEUE_PROFILING_ENABLE, NULL );
??? //創(chuàng)建三個OpenCL內(nèi)存對象，并把buf1的內(nèi)容通過隱式拷貝的方式
??? //buf1內(nèi)容拷貝到clbuf1,buf2的內(nèi)容通過顯示拷貝的方式拷貝到clbuf2
??? cl_mem clbuf1 = clCreateBuffer(context,
??????? CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
??????? BUFSIZE*sizeof(cl_float),buf1,
??????? NULL );

??? cl_mem clbuf2 = clCreateBuffer(context,
??????? CL_MEM_READ_ONLY ,
??????? BUFSIZE*sizeof(cl_float),NULL,
??????? NULL );

?? cl_event writeEvt;

??? status = clEnqueueWriteBuffer(queue, clbuf2, 1,
??????? 0, BUFSIZE*sizeof(cl_float), buf2, 0, 0, 0);

??? 上面這行代碼把buf2中的內(nèi)容拷貝到clbuf2,因為buf2位于host端，clbuf2位于device端，所以這個函數(shù)會執(zhí)行一次host到device的傳輸操作，或者說一次system memory到video memory的拷貝操作，所以我在該函數(shù)的后面放置了clFush函數(shù)，表示把command queue中的所有命令提交到device(注意：該命令并不保證命令執(zhí)行完成),所以我們調(diào)用函數(shù)waitForEventAndRelease來等待write緩沖的完成，waitForEventAndReleae是一個用戶定義的函數(shù)，它的內(nèi)容如下，主要代碼就是通過event來查詢我們的操作是否完成，沒完成的話，程序就一直block在這行代碼處，另外我們也可以用opencl中內(nèi)置的函數(shù)clWaitForEvents來代替clFlush和waitForEventAndReleae。

???? status = clFlush(queue);
???? //等待數(shù)據(jù)傳輸完成再繼續(xù)往下執(zhí)行
???? waitForEventAndRelease(&writeEvt);

??? cl_mem buffer = clCreateBuffer( context,
??????? CL_MEM_WRITE_ONLY,
??????? BUFSIZE * sizeof(cl_float),
??????? NULL, NULL );

????? kernel文件中放的是gpu中執(zhí)行的代碼，它被放在一個單獨的文件add.cl中，本程序中kernel代碼非常簡單，只是執(zhí)行兩個數(shù)組相加。kernel的代碼為：

?? //kernel文件為add.cl
??? const char * filename? = "add.cl";
??? std::string? sourceStr;
??? status = convertToString(filename, sourceStr);

convertToString也是用戶定義的函數(shù)，該函數(shù)把kernel源文件讀入到一個string中，它的代碼如下：

??? //創(chuàng)建程序?qū)ο?br /> ??? cl_program program = clCreateProgramWithSource(
??????? context,
??????? 1,
??????? &source,
??????? sourceSize,
??????? NULL);
??? //編譯程序?qū)ο?/span>
??? status = clBuildProgram( program, 1, &device, NULL, NULL, NULL );
??? if(status != 0)
??????? {
??????? printf("clBuild failed:%d\n", status);
??????? char tbuf[0x10000];
??????? clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, 0x10000, tbuf, NULL);
??????? printf("\n%s\n", tbuf);
??????? return -1;
??????? }

??? //創(chuàng)建Kernel對象
??? cl_kernel kernel = clCreateKernel( program, "vecadd", NULL );
??? //設(shè)置Kernel參數(shù)
??? cl_int clnum = BUFSIZE;
??? clSetKernelArg(kernel, 0, sizeof(cl_mem), (void*) &clbuf1);
??? clSetKernelArg(kernel, 1, sizeof(cl_mem), (void*) &clbuf2);
??? clSetKernelArg(kernel, 2, sizeof(cl_mem), (void*) &buffer);

注意：在執(zhí)行kernel時候，我們只設(shè)置了global work items數(shù)量，沒有設(shè)置group size，這時候，系統(tǒng)會使用默認(rèn)的work group size，通常可能是256之類的。

??? //執(zhí)行kernel,Range用1維，work itmes size為BUFSIZE
??? cl_event ev;
??? size_t global_work_size = BUFSIZE;
??? clEnqueueNDRangeKernel( queue,
??????? kernel,
??????? 1,
??????? NULL,
??????? &global_work_size,
??????? NULL, 0, NULL, &ev);
?? status = clFlush( queue );
?? waitForEventAndRelease(&ev);

??? //數(shù)據(jù)拷回host內(nèi)存
??? cl_float *ptr;

??? cl_event mapevt;
??? ptr = (cl_float *) clEnqueueMapBuffer( queue,
??????? buffer,
??????? CL_TRUE,
??????? CL_MAP_READ,
??????? 0,
??????? BUFSIZE * sizeof(cl_float),
??????? 0, NULL, NULL, NULL );

?? status = clFlush( queue );
?? waitForEventAndRelease(&mapevt);

???
??? //結(jié)果驗證，和cpu計算的結(jié)果比較
??? if(!memcmp(buf, ptr, BUFSIZE))
??????? printf("Verify passed\n");
??? else printf("verify failed");

??? if(buf)
??????? free(buf);
??? if(buf1)
??????? free(buf1);
??? if(buf2)
??????? free(buf2);

????? 程序結(jié)束后，這些opencl對象一般會自動釋放，但是為了程序完整，養(yǎng)成一個好習(xí)慣，這兒我加上了手動釋放opencl對象的代碼。

??? //刪除OpenCL資源對象
??? clReleaseMemObject(clbuf1);
??? clReleaseMemObject(clbuf2);
??? clReleaseMemObject(buffer);
??? clReleaseProgram(program);
??? clReleaseCommandQueue(queue);
??? clReleaseContext(context);
??? return 0;
??? }

程序執(zhí)行后的界面如下：

完整的代碼請參考：

工程文件gclTutorial1

代碼下載：

http://files.cnblogs.com/mikewolf2002/gclTutorial.zip

???? 首先我們建立工程文件gclTutorial2,在其中增加類gclFile,該類主要用來讀取文本kernel文件，或者讀寫二進(jìn)制kernel文件。

class gclFile
{
public:
??? gclFile(void);
??? ~gclFile(void);

??? //打開opencl kernel源文件(文本模式)
??? bool open(const char* fileName);

??? //讀寫二進(jìn)制kernel文件
??? bool writeBinaryToFile(const char* fileName, const char* birary, size_t numBytes);
??? bool readBinaryFromFile(const char* fileName);

…

}

gclFile中三個讀寫kernel文件的函數(shù)代碼為：

現(xiàn)在，在main.cpp中，我們就可以用gclFile類的open函數(shù)來讀入kernel源文件了：

//kernel文件為add.cl

gclFile kernelFile;
if(!kernelFile.open("add.cl"))
??? {
??? printf("Failed to load kernel file \n");
??? exit(0);
??? }
const char * source = kernelFile.source().c_str();
size_t sourceSize[] = {strlen(source)};
//創(chuàng)建程序?qū)ο?/span>
cl_program program = clCreateProgramWithSource(
??? context,
??? 1,
??? &source,
??? sourceSize,
??? NULL);

??? 編譯好kernel后，我們可以通過下面的代碼，把編譯好的kernel存儲在一個二進(jìn)制文件addvec.bin中，在教程4種，我們將會直接裝入這個二進(jìn)制的kernel文件。

??? 我們還會建立一個計時器類gclTimer,用來統(tǒng)計時間，這個類主要用QueryPerformanceFrequency得到時鐘頻率，用QueryPerformanceCounter得到流逝的ticks數(shù)，最終得到流逝的時間。函數(shù)非常簡單，

class gclTimer
{
public:
??? gclTimer(void);
??? ~gclTimer(void);

private:

??? double _freq;
??? double _clocks;
??? double _start;
public:
??? void Start(void); // 啟動計時器
??? void Stop(void); //停止計時器
??? void Reset(void); //復(fù)位計時器
??? double GetElapsedTime(void); //計算流逝的時間
};

下面我們在cpu端執(zhí)行數(shù)組加法時，增加計時器的代碼：

gclTimer clTimer;
clTimer.Reset();
clTimer.Start();
//cpu計算buf1,buf2的和
for(i = 0; i < BUFSIZE; i++)
??? buf[i] = buf1[i] + buf2[i];
clTimer.Stop();
printf("cpu costs time:%.6f ms \n ", clTimer.GetElapsedTime()*1000 );

同理在gpu執(zhí)行kernel代碼，以及copy gpu結(jié)果到cpu時候，增加計時器代碼：

最終程序執(zhí)行界面如下，在bufsize為262144時，在我的顯卡上gpu還有cpu快呢…，在程序目錄，我們可以看到也產(chǎn)生了vecadd.bin文件了。

完整的代碼請參考：

工程文件gclTutorial2

代碼下載：

http://files.cnblogs.com/mikewolf2002/gclTutorial.zip

本教程中，我們使用上一篇教程中產(chǎn)生的二進(jìn)制kernel文件vecadd.bin作為輸入來創(chuàng)建程序?qū)ο?#xff0c;程序代碼如下：

//kernel文件為vecadd.bin
gclFile kernelFile;
if(!kernelFile.readBinaryFromFile("vecadd.bin"))
??? {
??? printf("Failed to load binary file \n");
??? exit(0);
??? }
const char * binary = kernelFile.source().c_str();
size_t binarySize = kernelFile.source().size();

cl_program program = clCreateProgramWithBinary(context,
??? 1,
??? &device,
??? (const size_t *)&binarySize,
??? (const unsigned char**)&binary,
??? NULL,
??? NULL);

程序執(zhí)行的界面和教程3中一摸一樣…

完整的代碼請參考：

工程文件gclTutorial3

代碼下載：

http://files.cnblogs.com/mikewolf2002/gclTutorial.zip

//執(zhí)行kernel,Range用2維，work itmes size為width*height,
cl_event ev;
size_t globalThreads[] = {width, height};
size_t localx, localy;
if(width/8 > 4)
??? localx = 16;
else if(width < 8)
??? localx = width;
else localx = 8;

if(height/8 > 4)
??? localy = 16;
else if (height < 8)
??? localy = height;
else localy = 8;

size_t localThreads[] = {localx, localy};// localx*localy應(yīng)該是64的倍數(shù)
printf("global_work_size =(%d,%d), local_work_size=(%d, %d)\n",width,height,localx,localy);

clTimer.Reset();
clTimer.Start();
clEnqueueNDRangeKernel( queue,
??? kernel,
??? 2,
??? NULL,
??? globalThreads,
??? localThreads, 0, NULL, &ev);

注意：在上面代碼中，定義global threads以及l(fā)ocal threads數(shù)量，都是通過二維數(shù)組的方式進(jìn)行的。

??? 新的Kernel代碼如下：

????? 我們在kernel中增加了#pragma OPENCL EXTENSION cl_amd_printf : enable，以便在kernel中通過printf函數(shù)進(jìn)行debug，這是AMD的一個擴(kuò)展。printf還可以直接打印出float4這樣的向量，比如printf(“%v4f”, vec)。

????? 另外，在main.cpp中增加一行代碼：

//告訴driver dump il和isa文件
_putenv("GPU_DUMP_DEVICE_KERNEL=3");
????? 我們可以在程序目錄dump出il和isa形式的kernel文件，對于熟悉isa匯編的人，這是一個很好的調(diào)試performance的方法。

???? 在最新的app sdk 2.7中，在kernel中使用printf的時候，這個程序會hang在哪兒，以前沒這種情況。

程序執(zhí)行界面。

完整的代碼請參考：

工程文件gclTutorial4

代碼下載：

http://files.cnblogs.com/mikewolf2002/gclTutorial.zip

????? 首先我們建立一個gFreeImage類，用來裝入圖像，該類主要調(diào)用FreeImage的函數(shù)，首先會初始化FreeImage庫，然后根據(jù)文件名猜測圖像文件格式，最終load圖像文件到變量FIBITMAP *bitmap中去。同時，我們還定義了2個緩沖

unsigned char *imageData;
unsigned char *imageData4;

用來存放圖像數(shù)據(jù)，之所以定義imageData4,是因為通常的圖片文件，比如jpg,bmp都是3個通道，沒有包括alpha通道，但是在gpu中處理數(shù)據(jù)時候，通常以vector4或者vector的形式進(jìn)行，不能以vector3進(jìn)行，所以我們裝入圖像后，會把imageData指向圖像數(shù)據(jù)，同時生成包括alpha通道的圖像數(shù)據(jù)imageData4。

???? 另外，我們還定義了一個函數(shù)LoadImageGrey，該函數(shù)用來裝入灰度圖，灰度圖一個像素用一個uchar表示。

在main.cpp中，我們首先定義一個cpu處理圖像旋轉(zhuǎn)的函數(shù)：

//CPU旋轉(zhuǎn)圖像
void cpu_rotate(unsigned char* inbuf, unsigned char* outbuf, int w, int h,float sinTheta, float cosTheta)
??? {
??? int i, j;
??? int xc = w/2;
??? int yc = h/2;

??? for(i = 0; i < h; i++)
??????? {
??????? for(j=0; j< w; j++)
??????????? {
??????????? int xpos =? ( j-xc)*cosTheta - (i-yc)*sinTheta+xc;???
??????????? int ypos =? (j-xc)*sinTheta + ( i-yc)*cosTheta+yc;

??????????? if(xpos>=0&&ypos>=0&&xpos<w&&ypos<h)
??????????????? outbuf[ypos*w + xpos] = inbuf[i*w+j];
??????????? }
??????? }
??? }

??? 在main函數(shù)中，我們首先會裝入圖像文件，代碼如下：

??? 之后，定義2個cl memory對象，一個用來放原始圖像，一個用來放旋轉(zhuǎn)后的圖像。

//創(chuàng)建2個OpenCL內(nèi)存對象
cl_mem d_ip = clCreateBuffer(
??? context, CL_MEM_READ_ONLY,
??? mem_size,
??? NULL, NULL);
cl_mem d_op = clCreateBuffer(
??? context, CL_MEM_WRITE_ONLY,
??? mem_size,
??? NULL, NULL);

cl_event writeEvt;
status = clEnqueueWriteBuffer (???
??? queue , d_ip, CL_TRUE,
??? 0, mem_size, (void *)src_image,
??? 0, NULL, &writeEvt);
//等待數(shù)據(jù)傳輸完成再繼續(xù)往下執(zhí)行
status = clFlush(queue);
waitForEventAndRelease(&writeEvt);
//clWaitForEvents(1, &writeEvt);

?? 旋轉(zhuǎn)kernel函數(shù)需要傳入6個參數(shù)：

//創(chuàng)建Kernel對象
cl_kernel kernel = clCreateKernel( program, "image_rotate", NULL );
//設(shè)置Kernel參數(shù)
float sintheta = 1, costheta = 0;
clSetKernelArg(kernel, 0, sizeof(cl_mem),? (void *)&d_ip);
clSetKernelArg(kernel, 1, sizeof(cl_mem),? (void *)&d_op);
clSetKernelArg(kernel, 2, sizeof(cl_int),? (void *)&W);
clSetKernelArg(kernel, 3, sizeof(cl_int),? (void *)&H);
clSetKernelArg(kernel, 4, sizeof(cl_float), (void *)&sintheta);
clSetKernelArg(kernel, 5, sizeof(cl_float), (void *)&costheta);
kernel執(zhí)行的代碼為：

//執(zhí)行kernel,Range用2維，work itmes size為W*H,
cl_event ev;
size_t globalThreads[] = {W, H};
size_t localThreads[] = {16, 16}; // localx*localy應(yīng)該是64的倍數(shù)
printf("global_work_size =(%d,%d), local_work_size=(16, 16)\n",W,H);

clTimer.Reset();
clTimer.Start();
clEnqueueNDRangeKernel( queue,
??? kernel,
??? 2,
??? NULL,
??? globalThreads,
??? localThreads, 0, NULL, &ev);
//沒有設(shè)置local group size時候，系統(tǒng)將會自動設(shè)置為 (256,1)
status = clFlush( queue );
waitForEventAndRelease(&ev);
//clWaitForEvents(1, &ev);

clTimer.Stop();
printf("kernal total time:%.6f ms \n ", clTimer.GetElapsedTime()*1000 );

kernel函數(shù)代碼為：

gpu執(zhí)行完畢后，旋轉(zhuǎn)后的圖像保存在lenna_rotate.jpg,我們還會用cpu rotate函數(shù)執(zhí)行一次旋轉(zhuǎn)，同時把生成的圖像保存到cpu_lenna_rotate.jpg。

完整的代碼請參考：

工程文件gclTutorial5

代碼下載：

http://files.cnblogs.com/mikewolf2002/gclTutorial.zip

???? 灰度圖的histogram計算，首先要選擇bin(中文可以稱作槽)的數(shù)量，對于灰度圖，像素的范圍通常是[0-255],所以bin的數(shù)目就是256，然后我們循環(huán)整幅圖像，統(tǒng)計出每種像素值出現(xiàn)的次數(shù)，放到對應(yīng)的bin中。比如bin[0]中放的就是整幅圖像中灰度值為0的像素個數(shù)，bin[1]中放的就是整幅圖像中灰度值為1的像素個數(shù)……

???? 下面的直方圖就是灰度圖lenna對應(yīng)的直方圖。

????? 灰度圖直方圖的cpu計算特別簡單，定義一個數(shù)組hostBin[256],初始化所有數(shù)組元素為0，然后循環(huán)整幅圖像，得到直方圖，代碼如下：

//cpu求直方圖
void cpu_histgo()
??? {
??? int i, j;
??? for(i = 0; i < height; ++i)
??????? {
??????? for(j = 0; j < width; ++j)
??????????? {
??????????? //printf("data: %d\n", data[i * width + j] );
??????????? hostBin[data[i * width + j]]++;
??????????? //printf("hostbin %d=%d\n", data[i * width + j], hostBin[data[i * width + j]]);
??????????? }
??????? }
??? }

??? 如何使用opencl,來計算灰度圖，就沒有那么簡單了。我們知道gpu的優(yōu)勢是并行計算，如何把圖像分塊，來并行計算直方圖，是我們討論的重點。下面是一副512*512圖像的thread，workgroup劃分：

???? 我們設(shè)定圖像的寬度是bins的整數(shù)倍，即256的倍數(shù)，高度是workgroup size(本程序中，設(shè)置為128)的倍數(shù)，如果圖像高寬不是bins和workgroup size的倍數(shù)，則我們通過下面的公式把圖像的寬度和高度變成它們的倍數(shù)：

//width是binSize的整數(shù)倍，height是groupsize的整數(shù)倍
width = (width / binSize ? width / binSize: 1) * binSize;
height = (height / groupSize ? height / groupSize: 1) * groupSize;

???? 則512*512的圖像可以分為8個work group，每個workgroup包括128個thread，每個thread計算256個像素的直方圖，并把結(jié)果放到該thread對應(yīng)的local memroy空間，在kenrel代碼結(jié)束前，合并一個workgroup中所有thread的直方圖，生成一個workgroup塊的直方圖，最后在host端，合并8個workgroup塊的直方圖，產(chǎn)生最終的直方圖。

??? openCL的memory對象主要有3個，dataBuffer用來傳入圖像數(shù)據(jù)，而minDeviceBinBuf大小是workgroup number *256, 即每個workgroup對應(yīng)一個bin，另外一個kernel函數(shù)的第二個參數(shù)，它的大小為workgroup size*256, 用于workgroup中的每個線程存放自己256個像素的直方圖結(jié)果。

//創(chuàng)建2個OpenCL內(nèi)存對象
dataBuf = clCreateBuffer(
??? context,
??? CL_MEM_READ_ONLY,
??? sizeof(cl_uchar) * width? * height,
??? NULL,
??? 0);

//該對象存放每個block塊的直方圖結(jié)果
midDeviceBinBuf = clCreateBuffer(
??? context,
??? CL_MEM_WRITE_ONLY,
??? sizeof(cl_uint) * binSize * subHistgCnt,
??? NULL,
??? 0);

?? …

??? status = clSetKernelArg(kernel, 1, groupSize * binSize * sizeof(cl_uchar), NULL);//local memroy size, lds for amd

下面看看kernel代碼是如何計算workgroup塊的直方圖。

__kernel
void histogram256(__global const uchar* data,
????????????????? __local uchar* sharedArray,
????????????????? __global uint* binResult)
{
??? size_t localId = get_local_id(0);
??? size_t globalId = get_global_id(0);
??? size_t groupId = get_group_id(0);
??? size_t groupSize = get_local_size(0);

下面這部分代碼初始化每個thread對應(yīng)的local memory，也就是對應(yīng)的256個bin中計數(shù)清零。sharedArray大小是workgroup size * 256 = 128 * 256

??? //初始化共享內(nèi)存
??? for(int i = 0; i < BIN_SIZE; ++i)
??????? sharedArray[localId * BIN_SIZE + i] = 0;

通過barrier設(shè)置workgroup中所有thread的同步點，保證所有thread都完成初始化操作。

??? barrier(CLK_LOCAL_MEM_FENCE);
???

下面的代碼，計算thread中，256個像素的直方圖，比如對于workgroup 0中的thread 0它計算的256個像素為綠條的部分像素，注意：每個thread的包含的像素并不是連續(xù)的。

??? //計算thread直方圖
??? for(int i = 0; i < BIN_SIZE; ++i)
??? {
??????? uint value = (uint)data[groupId * groupSize * BIN_SIZE + i * groupSize + localId];
??????? sharedArray[localId * BIN_SIZE + value]++;
??? }
??? 通過fence，保證每個thread都完成各自的直方圖計算。
??? barrier(CLK_LOCAL_MEM_FENCE);?
??? 下面是合并各個thread的直方圖形成整個workgroup像素塊的直方圖，每個thread合并2個bin，比如thread 0，合并bin0和bin128。

?? //合并workgroup中所有線程的直方圖，產(chǎn)生workgroup直方圖
??? for(int i = 0; i < BIN_SIZE / groupSize; ++i)
??? {
??????? uint binCount = 0;
??????? for(int j = 0; j < groupSize; ++j)
??????????? binCount += sharedArray[j * BIN_SIZE + i * groupSize + localId];
???????????
??????? binResult[groupId * BIN_SIZE + i * groupSize + localId] = binCount;
??? }
}

最終在host端，我們還要把每個workgroup塊的直方圖合并成得到整個圖像的直方圖，主要代碼如下：

// 合并子塊直方圖值

for(i = 0; i < subHistgCnt; ++i)
??? {
??? for( j = 0; j < binSize; ++j)
??????? {
??????? deviceBin[j] += midDeviceBin[i * binSize + j];
??????? }
??? }

完整的代碼請參考：

工程文件gclTutorial7

代碼下載：

http://files.cnblogs.com/mikewolf2002/gclTutorial.zip

現(xiàn)在我們利用上一篇教程的方法，來統(tǒng)計一副RGBA圖像中，有多少個像素點（該像素點滿足R, G, B, A任意分量>=5)？我考慮的方法是建立256 bin的直方圖，對于一個像素，求max(R, G,B,A),用該值決定該像素點進(jìn)入那個bin，這樣求出直方圖后，width*height - hostBin[0] - hostBin[1] - hostBin[2] - hostBin[3] - hostBin[4]，即為我們要的結(jié)果。

???? 本教程代碼基本上和上一篇教程中代碼是一樣的，區(qū)別主要包括以下2點：

1、我們裝入的是RGBA 彩色圖。

//裝入圖像
unsigned char *src_image=0;
gFreeImage img;
if(!img.LoadImage("../lenna.jpg"))
??? {
??? printf("can‘t load lenna.jpg\n");
??? exit(0);
??? }
src_image = img.getImageData(width,height);

2、在kernel代碼中，有點小變化，kernel代碼如下：

?

完整的代碼請參考：

工程文件gclTutorial8

代碼下載：

http://files.cnblogs.com/mikewolf2002/gclTutorial.zip

??? 首先我們用amd的opencl profiler分析一下程序性能(不會找不到用吧，點擊view-other windows-app profiler…,然后就看到了…)。

???? 下面我們來分析我們的kernel代碼中memory操作：

???? 首先是shared memory的的初始化，我們知道shared memory是local memory，被一個workgroup中的所有thread，或者說work item共享。在amd硬件系統(tǒng)中，local memory是LDS，它通常是為32k，分為32個bank，dword字節(jié)地址，每個bank 512個item，我們可以通過函數(shù)得到自己系統(tǒng)中的local memory數(shù)量：

cl_ulong DeviceLocalMemSize;
clGetDeviceInfo(device,
??? CL_DEVICE_LOCAL_MEM_SIZE,
??? sizeof(cl_ulong),
??? &DeviceLocalMemSize,
??? NULL);

lds的示意圖如下，對于每個bank，同時只能有一個讀寫請求，如果兩個thread都讀寫bank1，那個必須串行訪問，這就稱作bank conflict。

kernel初始化local memory的代碼如下：

//初始化共享內(nèi)存
for(int i = 0; i < BIN_SIZE; ++i)
??? sharedArray[localId * BIN_SIZE + i] = 0;

???? 在同一時間，thread0訪問地址0(bank1),thread1,訪問地址256，也在bank1，…,這樣就有很多bank conflit，降低程序的性能。從profiler里面可以看到，lds bank conflit為13.98，很高的比例，所以此時同時運行的thread就比較少，只有總wave(每個wave 64個thread)的12%（我曾經(jīng)默認(rèn)lds內(nèi)存分配是0，這樣我們就省去了這些代碼，但是實際上分配內(nèi)存是一些隨機(jī)的值…)。

第二段memory操作的代碼為：

??? //計算thread直方圖
??? for(int i = 0; i < BIN_SIZE; ++i)
??? {
??????? uint value = (uint)data[groupId * groupSize * BIN_SIZE + i * groupSize + localId];
??????? sharedArray[localId * BIN_SIZE + value]++;
??? }

其中有l(wèi)ds的操作，也有g(shù)lobal memory的操作，對于全局memory的訪問，在同一時刻，thread0訪問i=0的memory

，thread1訪問相鄰的memory單元…,這是對于global memory的訪問會采用合并讀寫的方式(coalencing)，就是一個memory請求返回16個dword，也就是一個請求滿足16個thread，提高memory利用率。此時對lds的寫是隨機(jī)的，根據(jù)value的值決定，不能控制…

?

最后一段memory讀寫的代碼：

//合并workgroup中所有線程的直方圖，產(chǎn)生workgroup直方圖
for(int i = 0; i < BIN_SIZE / groupSize; ++i)
{
???? uint binCount = 0;
???? for(int j = 0; j < groupSize; ++j)
???????? binCount += sharedArray[j * BIN_SIZE + i * groupSize + localId];
????????
???? binResult[groupId * BIN_SIZE + i * groupSize + localId] = binCount;
}

其中l(wèi)ds的讀寫如下圖，此時每個線程訪問不同的bank，因為amd lds訪問就是以32為單位，所以實際上，這段代碼不會有bank conflit。

???? 我們的矩陣是一副二維灰度圖像256*256，矩陣的轉(zhuǎn)置也就是圖像的轉(zhuǎn)置。每個thread處理16(4*4)個pixel(uchar)，workgroup的size是(16,16)。

???? 下面直接看shader代碼：

uint wiWidth? = get_global_size(0);

uint gix_t = get_group_id(0);
uint giy_t = get_group_id(1);???
uint num_of_blocks_x = get_num_groups(0);

uint giy = gix_t;
uint gix = (gix_t+giy_t)%num_of_blocks_x;

uint lix = get_local_id(0);
uint liy = get_local_id(1);

uint blockSize = get_local_size(0);

uint ix = gix*blockSize + lix;
uint iy = giy*blockSize + liy;
int index_in = ix + (iy)*wiWidth*4;

// 通過合并讀寫把輸入數(shù)據(jù)裝入到lds中

int ind = liy*blockSize*4+lix;

block[ind]??????? = input[index_in];
block[ind+blockSize]??? = input[index_in+wiWidth];
block[ind+blockSize*2] = input[index_in+wiWidth*2];
block[ind+blockSize*3] = input[index_in+wiWidth*3];

???? 因為workgroup size是（16,16），所以lix,liy的取值范圍都是0-15，下面我們通過圖片看下，lix=0 liy=0,lix=1 liy=0時候，ind,以及index_in的值，從而得到輸入圖像數(shù)據(jù)如何映射到local memory中。

lix=0 liy=0

lix=1 liy=0

????? 下面是影射關(guān)系，(0,0) thread處理的16個pixel用血紅色表示，它們映射到lds的0bank和16bank，(1,0)thread處理的像素用綠色表示，它們映射到lds的bank1和bank17，有效的避免了bank conflit，而全局memory的訪問不同thread對應(yīng)連續(xù)的全局memory空間，可以實現(xiàn)合并讀寫，從而提高程序性能。

? 把轉(zhuǎn)置的數(shù)據(jù)寫到全局memory中的代碼如下：

ix = giy*blockSize + lix;
iy = gix*blockSize + liy;
int index_out = ix + (iy)*wiWidth*4;

ind = lix*blockSize*4+liy;
uchar4 v0 = block[ind];
uchar4 v1 = block[ind+blockSize];
uchar4 v2 = block[ind+blockSize*2];
uchar4 v3 = block[ind+blockSize*3];

// 通過合并讀寫把lds中數(shù)據(jù)寫回到全局memory中
output[index_out]??????????? = (uchar4)(v0.x, v1.x, v2.x, v3.x);
output[index_out+wiWidth]??? = (uchar4)(v0.y, v1.y, v2.y, v3.y);
output[index_out+wiWidth*2]??? = (uchar4)(v0.z, v1.z, v2.z, v3.z);
output[index_out+wiWidth*3]??? = (uchar4)(v0.w, v1.w, v2.w, v3.w);
對應(yīng)copy關(guān)系圖如下：

完整的代碼請參考：

工程文件gclTutorial9

代碼下載：

稍后提供

???? 在程序中，我們設(shè)置workgroup size為256，kernel的輸入、輸出緩沖參數(shù)都用uint4的格式，這樣我們原始求和的數(shù)組大小為256*4的倍數(shù)，數(shù)據(jù)類型為uint。我們設(shè)定每個workgroup處理處理512個uint，即2048個uint

???? 為了簡便期間，我們輸出數(shù)組長度定為4096，即需要2個workgruop來處理。

???

kernel代碼如下：

__kernel void reduce(__global uint4* input, __global uint4* output, __local uint4* sdata)
{
??? // 把數(shù)據(jù)裝入lds
??? unsigned int tid = get_local_id(0);
??? unsigned int bid = get_group_id(0);
??? unsigned int gid = get_global_id(0);

??? unsigned int localSize = get_local_size(0);
??? unsigned int stride = gid * 2;
??? sdata[tid] = input[stride] + input[stride + 1];

??? barrier(CLK_LOCAL_MEM_FENCE);
??? // 在lds中進(jìn)行reduction操作，得到數(shù)組求和的結(jié)果
??? for(unsigned int s = localSize >> 1; s > 0; s >>= 1)
??? {
??????? if(tid < s)
??????? {
??????????? sdata[tid] += sdata[tid + s];
??????? }
??????? barrier(CLK_LOCAL_MEM_FENCE);
??? }

?? // 把一個workgroup計算的結(jié)果輸出到輸出緩沖，是一個uint4，還需要在host端再進(jìn)行一次reduction過程
??? if(tid == 0) output[bid] = sdata[0];
}

???? 在程序中，global和local的NDRange，我們都用一維的形式。下面以圖的方式看下kernel代碼是如何執(zhí)行的：

??? ? 對第一個workgroup中的第一個thread的來說，它首先進(jìn)行一次reduction操作，把兩個uint4相加，放到lds(shared memory)中，然后再在lds中進(jìn)行reduction操作，此時要從global memory中取數(shù)據(jù)，可以看出連續(xù)的thread訪問連續(xù)的global memory，這時可以利用合并讀寫。

????? 申請的shared memory大小為groupsize*sizeof(uint4)，相加后uint4放入32bank的lds中，放置的方式應(yīng)該是如下圖所示，因為放入的是uint4，所以會放入連續(xù)的4個bank中(每個bank都是dword寬)，可見只能同時有8個thread訪問lds，所以會有一定程序的bank conflit。從App profiler session，我們可以看到：

????? 接下來，kernel會通過一個for循環(huán)迭代執(zhí)行reduction操作，求得一個workgroup中的uint4的和。

迭代的第一次s=128,這時會執(zhí)行如下圖的兩兩相加，workgroup中同時執(zhí)行的thread為128,thread local id大于等于128的線程都不會做什么事情，在每個循環(huán)的末尾，有一個barrier來同步所有thread，以便所有thread都完成這次循環(huán)后再進(jìn)入下一次循環(huán)。

? ?? 第二次迭代的時候，只剩下前面128個uint4，workgroup中同時執(zhí)行的thread為64。最后，當(dāng)s=1時候，完成迭代reduction操作，然后把thread0（第一個thread)的結(jié)果輸出。

???? 在host段，我們還要做一次相加操作，把不同workgroup得到的uint4，拆分成uint，并相加求得最終的結(jié)果。

//在cpu reduction各個workgroup的結(jié)果以及uint4分量 reduction
output = 0;
for(int i = 0; i < numBlocks * VECTOR_SIZE; ++i)
??? output += outMapPtr[i];
printf("gpu reduction result:%d\n", output);
if(refOutput==output) printf("passed\n");

程序執(zhí)行后結(jié)果如下：

完整的代碼請參考：

工程文件gclTutorial11

代碼下載：

稍后提供

?

總結(jié)

以上是生活随笔為你收集整理的AMD OpenCL 大学课程的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。