1. 程式人生 > >CUDA GPU高效能運算

CUDA GPU高效能運算

分享一下我老師大神的人工智慧教程!零基礎,通俗易懂!http://blog.csdn.net/jiangjunshow

也歡迎大家轉載本篇文章。分享知識,造福人民,實現我們中華民族偉大復興!

               

序言

CUDA是異構程式設計的一個大頭,洋洋灑灑的看了些資料,但是,感覺這個技術沒有像C++或者Java那樣有自己的權威的《程式設計思想》來指導系統學習,總是感覺心裡不踏實,是不是自己還沒掌握深入、或者說心裡沒底氣說自己已經入門了、已經熟悉了、已經精通了。站在一個初學者的角度,作為一個筆記式的記錄,講解自己學習和理解

CUDA過程中的一些列想到的、碰到的問題。享受一個東西不一定是結果,可以是從無知到了解到精通的這個整個過程。

給自己提幾個問題

對的,我想要做什麼事情的時候,習慣性的給自己提如下問題:

問題1GPU高效能運算之CUDA(下午簡稱CUDA技術)是幹嘛用的,我為什麼要學它?

這個問題如果我回答是,純粹為了掌握一門新的技術,如果你還是學生則可以,如果你是一個手裡有專案的工程師,那麼,我覺得沒什麼必要去學這個東西。我個人理解CUDA是計算機裡面的邊緣技術,是對程式執行效能的提高的一種方式,可以理解成一個工具。工具這種東西,你需要用的時候再去學,你不需要用的時候,你知道有這回事情就可以了。如果你用不到的東西,你也很貪心什麼都去學,第一學不精,第二計算機的技術太多太廣,對應我這樣的智商一般的人來說是不靠譜的。

我要學習CUDA,因為專案的任務可分解性較大,粒度相關性小,某專案大概一個計算任務可以分解成3000*3000*300規模的計算,而且這些元素之間完全獨立。又加入,你需要對全市100萬人,每個人計算下該個體的每年的收入與支出的淨值(從1900年計算到2000年),淡然有人說放excel中計算不就得了,好吧……我只是講個通俗點的例子。那麼這個計算規模是1000000*10,倘若要對每個人每年在做點其他什麼高階的演算法得出一個什麼指數,恐怕excel還是實現起來不太容易。所以,這種並行度很大的問題,我們可以考慮用CUDA來解決。也許你也知道,搞影象的,CUDA就顯得很重要啦。


問題2:我的基礎是什麼?

學習新的技術,我喜歡和之前學過的某個還熟悉的東西做為比對,這樣理解起來可能會快一些。比如學Java

的時候,我想著C++,學UML的時候,我想著面向物件程式設計,學CUDA呢?我有什麼麼?可能很多人會沒有異構程式設計的歷史。我很幸運,之前弄過半年的OpenMP,對的——多核程式設計技術。後來我瞭解到,如果是多個GPU,是需要用到OpenMP來做的。OpenMP本身和CUDA好像並沒什麼關係,但是,片上多核的並行演算法是很想通的。很好的一個例子是:奇偶排序。哈哈,給自己提了點學好CUDA的信心。

好了,明白自己的需求和基礎,給自己一個學習的定位,什麼地方該花時間去琢磨,相比應該很清楚。

你好,GPU

GPU有兩大開發商——英偉達和AMD,支援CUDA的是英偉達,很好啊,之前買電腦是有先見之明的,買了英偉達的顯示卡——GT520,不是特別高階,和單位的GTX650ti2G比起來,有那麼一點遜色,但是,好歹可以跑CUDA

英偉達支援CUDA程式設計的顯示卡型號從G8800開始,都是可以的。一開始作為影象處理用,而今,天文地理、數學金融、醫療軍事等等,都開始嘗試發揮GPU的優勢。

GPU的計算核心也是隔年換代,現在已經倒GK10X了,計算能力也是逐漸提升,目前已經最好的有3.0GPU的架構從原始的到費米的,再到開普勒的,我們沒必要去一個個瞭解,我們先了解GPU的這個大概的歷史,免得和人交談說不出一和二,脫離菜鳥嘛!

GPUCPU的區別可以參考下,講的還算詳細和明瞭。

http://www.cnblogs.com/viviman/archive/2012/11/26/2789113.html

可以這樣的去理解,單核CPU多執行緒並行,是感官上的並行,世界上在CPU上還是序列指令在跑;而在GPU上,才叫真正的並行!需要介紹下CUDA1.0開發包是支援在CPU上模擬GPU開發的,其原理就是用多執行緒來模擬;而現在的版本就不支援了。最新的是5.0的,官網上是下不到之前的了,反正我是找到腿軟了還沒找到。1.0是古董了!如果你有,一定要給我開開眼見哦。

我們的CUDA程式設計,很明顯是GPUCPU一起來處理的嘛——異構程式設計!對的,我自己的理解是就一個工程而言:CPU處理序列計算業務,GPU處理平行計算業務,這裡將的並行都是並行度可觀的哦,不是說兩個元素你也來GPU上計算,那樣是不環保的——會浪費很多GPU的資源!

說道環保,我想多說一句,在GPU程式設計,很體現“綠色”理念。48個核,你如果寫的好,48個核全在幹活,而且乾的是有意義的活,那麼你是合理的利用資源,如果你只讓一個核在幹有意義的活,其他的都在空轉,那麼你很浪費電哦。

對於CPUGPU分配自己的業務,稍微畫個圖失意一下,如圖1

圖1 工程中GPUCPU的分工

總的來說,GPU只是幹計算並行度高的功能模組的活,一定不可以越權啦!

你好,CUDA

3.1 開發環境配置

第一次和同事交談,我說你和我說說C-U-D-A。他說:哭打……苦打…… ……。半天后,我說,你和我說說C-U-D-A,你說哭打是什麼東西……它說哭打就是C-U-D-A。我操,頓時傻眼了,原來這東西行業裡年哭打,對的,Cu - Da,連起來就是這麼發音的,我的無知啊。我們還是用中文解釋吧:CUDA的意思是統一計算裝置架構。

CUDA的整合開發環境可以參考下:

http://www.cnblogs.com/viviman/archive/2012/11/05/2775100.html

win7+vs2008+CUDA5.0的環境下體驗,是一種新的嘗試,我自己配的時候,很少或者就沒有5.0的配置部落格文件等等。因為5.0的那一場雪比2010年來的更晚一些……

5.0是和2.03.04.0都會有那麼一點不同的,是集成了SDKTOOL兩個東西,之前是分開的,現在是合在一起的。其實是差不多的,但是,這一分一合,就會給人很不習慣。不過,我雖然愚鈍,但是,試了幾次之後還是摸索成功。當第一個helloworld輸出後,心裡是有那麼一點小激動的,立馬跑到小區門口,買了半斤羊肉吃了,因為為了配這一套環境,我中午飯都沒吃!我這隻哭打小菜鳥就是可憐啊!

3.2 特殊的"hello world"

搭建好環境,你肯定想看看我的helloworld程式,對的,但是我不想輸出helloworld,我想幹一件事情是:我在CPU上建立個變數,傳到GPU中,然後,在GPU中賦值,然後傳出來。這件事情如果成功了,是不是可以說明,通了+GPU工作了!網上一搜的CUDAhelloworld程式,都是在CPU上輸出helloworld,那多不過癮。

__global__ void hello(char *ch){ch = {'h', 'e', 'l', 'l', 'o'};}int main(){……hello<<<1,1>>>(dev_ch);……return 0;}

I think you know my idea.

在你網上搜索到的程式的基礎上,作這樣的一個改變,相信自己動手的才是快樂的。

敲開程式設計的門

我習慣性的喜歡先看一門語言的關鍵字,CUDA的關鍵字很簡單很少:

函式型別

__global__: 用來修飾核心函式的,核心函式是什麼呢,核心函式是跑在GPU上的函式;與之對應的是主機函式,用__host__修飾,也可以預設,跑在CPU上。因此,CPU也叫主機,GPU也叫裝置。通常定義這個核心函式,我喜歡在函式名前加個kernel作為修飾,讓自己清楚點。

比如__global__ void kerneladd(float *a){}

__device__: 也是用來修飾核心函式的,那和__global__有什麼區別嗎?對的。__global__修飾的核心函式只能被主機函式呼叫;__device__修飾的核心函式只能被核心函式呼叫,應該很好理解。

__host__: 主機函式,供主機函式呼叫,可預設哦,一般情況下,都是預設的,知道這個東西就行。

儲存型別

暫存器:在核函式內 int i即表示暫存器變數。

__global__:     全域性記憶體。在主機函式中開闢和釋放。

__shared__:   共享儲存,每個block內的執行緒共享這個儲存。

__constant__:常量儲存,只讀。定義在所有函式之外,作用範圍整個檔案。

__texture__:   紋理儲存,只讀。記憶體不連續。

內建變數

dim3, threadId, blockId,  gridId

5 GPU也不允許偏心

並行的事情多了,我們作為GPU的指令分配者,不能偏心了——給甲做的事情多,而乙沒事做,個麼甲肯定不爽的來。所以,在GPU中,叫做執行緒網路的分配。首先還是來看下GPU的執行緒網路吧,圖2

執行緒網路

我們將具體點的,在主機函式中如果我們分配的是這樣的一個東西:

dim3 blocks(32, 32);

dim3 threads(16, 16);

dim3是神馬?dim3是一個內建的結構體,和linux下定義的執行緒結構體是個類似的意義的東西,dim3結構變數有xyz,表示3維的維度。不理解沒關係,慢慢看。

kernelfun<<<blocks, threads>>>();

我們呼叫kernelfun這個核心函式,將blocksthreads傳到<<<,>>>裡去,這句話可牛逼大了——相當於發號施令,命令那些執行緒去幹活。這裡使用了32*32 * 16*16個執行緒來幹活。你看明白了嗎?blocks表示用了二維的32*32block組,而每個block中又用了16*16的二維的thread組。好吧,我們這個施令動用了262144個執行緒!我們先不管GPU內部是如何排程這些執行緒的,反正我們這一句話就是用了這麼多執行緒。

那我們的核心函式kernelfun()如何知道自己執行的是哪個執行緒?這就是執行緒網路的特點啦,為什麼叫網路,是有講究的,網路就可以定格到網點:

比如int tid = threadId.x + blockId.x * 16

這裡有一個講究,block是有維度的,一維、二維、三維。

對於一維的blocktid = threadId.x

對於(DxDy)二維的blocktid = threadId.x + Dx*threadId.y

對於(DxDyDz)三維的blocktid = threadId.x + Dx*threadId.y + Dz*Dy*threadId.z

我習慣的用這樣的模式去分配,比較通用:

dim3 dimGrid();

dim3 dimBlock();

kerneladd<<<dimGrid, dimBlock>>>();

這可是萬金油啊,你需要做的事情是填充dimGriddimBlock的結構體建構函式變數,比如,dimGrid(16, 16)表示用了16*16的二維的block執行緒塊。

0,0)(0,1)(0,2)……(0,15

1,0)(1,1)(1,2)……(1,15

2,0)(2,1)(2,2)……(2,15

……

15,0)(15,1)(15,2)……(15,15

(,)是(dimGrid.x, dimGrid.y)的網格編號。

我們這麼理解吧,現在又一群人,我們分成16*16個小組(block),排列好,比如第3行第4列就指的是(2,3)這個小組。

dimBlock16,16)表示每個小組有16*16個成員,如果你想點名第3行第4列這個小組的裡面的第3行第4列那個同學,那麼,你就是在(2,3)這個block中選擇了(2,3)這個執行緒。這樣應該有那麼一點可以理解進去的意思了吧?不理解透徹麼什麼關係,這個東西本來就是cuda中最讓我糾結的事情。我們且不管如何分配執行緒,能達到最優化,我們的目標是先讓GPU正確地跑起來,計算出結果即可,管他高效不高效,管他環保不環保。

嘮叨了這麼多,下面我們用一個最能說明問題的例子來進一步理解執行緒網路分配機制來了解執行緒網路的使用。

一維網路執行緒

egint arr[1000],對每個陣列元素進行加1操作。

idea:我們最直接的想法,是排程1000個執行緒去幹這件事情。

first pro:我想用一個小組的1000個人員去幹活。這裡會存在這樣一個問題——一個小組是不是有這麼多人員呢?是的,這個事情你必須瞭解,連自己組內多少人都不知道,你也不配作指揮官呀。對的,這個引數叫做maxThreadsPerBlock,如何取得呢?

好吧,cuda定義了一個結構體cudaDeviceProp,裡面存入了一系列的結構體變數作為GPU的引數,出了maxThreadsPerBlock還有很多資訊哦,我們用到了再說。

maxThreadsPerBlock這個引數值是隨著GPU級別有遞增的,早起的顯示卡可能512個執行緒,我的GT520可以跑1024個執行緒,辦公室的GTX650ti2G可以跑1536個,無可非議,當然多多益善。一開始,我在想,是不是程式將每個block開的執行緒開滿是最好的呢?這個問題留在以後在說,一口吃不成胖子啦。

好吧,我們的陣列元素1000個,是可以在一個block中幹完的。

核心函式:

#define N 1000__gloabl__ void kerneladd(int *dev_arr)int tid = threadId.x; if (tid < 1000) dev_arr[tid] ++;}int main()int *arr, *dev_arr;// 習慣的我喜歡在核心函式引數變數前加個dev_作為標示 // 開闢主機記憶體,arr = (int*)malloc(N*sizeof(int)); // 開闢裝置記憶體 // 主機拷貝到裝置 kerneladd<<<1, N>>>(dev_arr); // 裝置拷貝到主機 // 列印 // 釋放裝置記憶體 // 釋放主機記憶體 return 0;}

呀,原來這麼簡單,個麼CUDA也忒簡單了哇!這中想法是好的,給自己提高信心,但是這種想法多了是不好的,因為後面的問題多了去了。

盆友說,1000個元素,還不如CPU來的快,對的,很多情況下,資料量並行度不是特別大的情況下,可能CPU來的更快一些,比較裝置與主機之間互相排程操作,是會有額外開銷的。有人就問了,一個10000個元素的陣列是不是上面提供的idea就解決不了啦?對,一個block人都沒怎麼多,如何完成!這個情況下有兩條路可以選擇——

第一,我就用一個組的1000人來幹活話,每個人讓他幹10個元素好了。

這個解決方案,我們需要修改的是核心函式:

__global__ void kernelarr(int *dev_arr)int tid = threadId.x; if(tid < 1000) // 只用0~999號執行緒 {  //每個執行緒處理10個元素,比如0號執行緒處理0、1001、2001、……9001  for(int i = tid; i<N; i=i+1000)  {  dev_arr[tid] ++;  } }}

第二,我多用幾個組來幹這件事情,比如我用10個組,每個組用1000人。

這個解決方案就稍微複雜了一點,注意只是一點點哦~因為,組內部怎麼幹活和最原始的做法是一樣的,不同之處是,我們調遣了10個組去幹這件事情。

首先我們來修改我們的主機函式:

int main(){……kerneladd<<<10, 1000>>>(dev_arr);//我們調遣了10個組,每個組用了1000人……}

盆友要問了,10個組每個組1000人,你怎麼點兵呢?很簡單啊,第1組第3個執行緒出列,第9組第9個執行緒出列。每個人用組號和組內的編號定了位置。線上程網路中,blockId.xthreadId.x就是對應的組號和組內編號啦,我必須要這裡開始形象點表示這個對應關係,如果這個對應關係是這樣子的[blockId.xthreadId.x],那麼我們的陣列arr[10000]可以這樣分配給這10個組去幹活:

(0,0)——arr[0](0,1)——arr[1],……(0,999)——arr[999]

(1,0)——arr[0+1*1000](1,1)——arr[1+1*1000],……(1,999)——arr[999+1*1000]

……

(9,0)——arr[0+9*1000](9,1)——arr[1+9*1000],……(9,999)——arr[999+9*1000]

是不是很有規律呢?對的,用blockId.xthreadId.x可以很好的知道哪個執行緒幹哪個元素,這個元素的下表就是threadId.x + 1000*blockId.x

這裡我想說的是,如果我們哪天糊塗了,畫一畫這個對應關係的表,也許,就更加清楚的知道我們分配的執行緒對應的處理那些東西啦。 

一維執行緒網路,就先學這麼多了。

二維網路執行緒

eg2int arr[32][16]二維的陣列自增1

第一個念頭,開個32*16個執行緒好了哇,萬事大吉!好吧。但是,朕現在想用二維執行緒網路來解決,因為朕覺得一個二維的網路去對映一個二維的陣列,朕看的更加明瞭,看不清楚自己的士兵,如何帶兵打仗!

我還是畫個對映關係:

一個block中,現在是一個二維的thread網路,如果我用了16*16個執行緒。

(0,0)(0,1),……(0,15)

(1,0)(1,1),……(1,15)

……

(15,0)(15,1),……(15,15)

呀,現在一個組內的人稱呼變了嘛,一維網路中,你走到一個小組裡,叫3號出列,就出來一個,你現在只是叫3號,沒人會出來!這個場景是這樣的,現在你班上有兩個人同名的人,你只叫名,他們不知道叫誰,你必須叫完整點,把他們的姓也叫出來。所以,二維網路中的(0,3)就是原來一維網路中的3,二維中的(i,j)就是一維中的(j+i*16)。不管怎麼樣,一個block裡面能處理的執行緒數量總和還是不變的。

一個grid中,block也可以是二維的,一個block中已經用了16*16thread了,那我們一共就32*16個元素,我們用2block就行了。

先給出一個程式碼清單吧,程式設計師都喜歡看程式碼,這段程式碼是我抄襲的。第一次這麼完整的放上程式碼,因為我覺得這個程式碼可以讓我說明我想說的幾個問題:

第一,二維陣列和二維指標的聯絡。

第二,二維執行緒網路。

第三,cuda的一些記憶體操作,和返回值的判斷。

#include <stdio.h> #include <stdlib.h> #include <cuda_runtime.h>  #define ROWS 32 #define COLS 16 #define CHECK(res) if(res!=cudaSuccess){exit(-1);} __global__ void Kerneltest(int **da, unsigned int rows, unsigned int cols) {  unsigned int row = blockDim.y * blockIdx.y + threadIdx.y;  unsigned int col = blockDim.x * blockIdx.x + threadIdx.x;  if (row < rows && col < cols)  {   da[row][col] = row * cols + col;  }  }  int main(int argc, char **argv) {  int **da = NULL;  int **ha = NULL;  int *dc = NULL;  int *hc = NULL;  cudaError_t res;  int r, c;  bool is_right = true;   res = cudaMalloc((void**)(&da), ROWS * sizeof(int*)); CHECK(res)  res = cudaMalloc((void**)(&dc), ROWS * COLS * sizeof(int)); CHECK(res)  ha = (int**)malloc(ROWS * sizeof(int*));  hc = (int*)malloc(ROWS * COLS * sizeof(int));   for (r = 0; r < ROWS; r++)  {   ha[r] = dc + r * COLS;  }  res = cudaMemcpy((void*)(da), (void*)(ha), ROWS * sizeof(int*), cudaMemcpyHostToDevice); CHECK(res)  dim3 dimBlock(16, 16);  dim3 dimGrid((COLS + dimBlock.x - 1) / (dimBlock.x), (ROWS + dimBlock.y - 1) / (dimBlock.y));  Kerneltest<<<dimGrid, dimBlock>>>(da, ROWS, COLS);  res = cudaMemcpy((void*)(hc), (void*)(dc), ROWS * COLS * sizeof(int), cudaMemcpyDeviceToHost); CHECK(res)   for (r = 0; r < ROWS; r++)  {   for (c = 0; c < COLS; c++)   {    printf("%4d ", hc[r * COLS + c]);    if (hc[r * COLS + c] != (r * COLS + c))    {     is_right = false;    }   }   printf("\n");  }  printf("the result is %s!\n", is_right ? "right" : "false");  cudaFree((void*)da);  cudaFree((void*)dc);  free(ha);  free(hc);  getchar();  return 0;  } 

簡要的來學習一下二維網路這個知識點, 

dim3 dimBlock(16, 16); 

//定義block內的thread二維網路為16*16

dim3 dimGrid((COLS + dimBlock.x - 1) / (dimBlock.x),  (ROWS + dimBlock.y - 1) / (dimBlock.y)); 

//定義grid內的block二維網路為1*2

unsigned int row = blockDim.y * blockIdx.y  +  threadIdx.y; 

//二維陣列中的行號

unsigned int col = blockDim.x * blockIdx.x  +  threadIdx.x; 

//二維執行緒中的列號

三維網路執行緒

dim3定義了三維的結構,但是,貌似二維之內就能處理很多事情啦,所以,我放棄學習三維。網上看到的不支援三維網路是什麼意思呢?先放一放。

給自己充充電

同一塊顯示卡,不管你是二維和三維或一維,其計算能力是固定的。比如一個block能處理1024個執行緒,那麼,一維和二維執行緒網路是不是處理的執行緒數一樣呢?

 

回答此問題,先給出網路配置的引數形式——<<<Dg,Db,Ns,S>>>,各個引數含義如下:

Dg:定義整個grid的維度,型別Dim3,但是實際上目前顯示卡支援兩個維度,所以,dim3<<Dg.x, Dg.y, 1>>>z維度預設只能為1,上面顯示出這個最大有65536*65536*1,每行有65536block,每列有65536block,整個grid中一共有65536*65536*1個block

Db:定義了每個block的維度,型別Dim3,比如512*512*64,這個可以定義3維尺寸,但是,這個地方是有講究了,三個維度的積是有上限的,對於計算能力1.01.1GPU,這個值不能大於768,對於1.21.3的不能大於1024,對於我們試一試的這塊級別高點的,不能大於1536。這個值可以獲取哦——maxThreadsPerBlock

Ns:這個是可選引數,設定最多能動態分配的共享記憶體大小,比如16k,單不需要是,這個值可以省略或寫0

S也是可選引數,表示流號,預設為0流這個概念我們這裡不說。

接著,我想解決幾個你肯定想問的兩個問題,因為我看很多人想我這樣的問這個問題:

1 block內的thread我們是都飽和使用嗎?

答:不要,一般來說,我們開128256個執行緒,二維的話就是16*16

2 grid內一般用幾個block呢?

答:牛人告訴我,一般來說是你的流處理器的4倍以上,這樣效率最高。

回答這兩個問題的解釋,我想抄襲牛人的一段解釋,解釋的好的東西就要推廣呀:

GPU的計算核心是以一定數量的Streaming Processor(SP)組成的處理器陣列,NV稱之為Texture Processing Clusters(TPC),每個TPC中又包含一定數量的Streaming Multi-Processor(SM),每個SM包含8SPSP的主要結構為一個ALU(邏輯運算單元),一個FPU(浮點運算單元)以及一個Register File(暫存器堆)SM內包含有一個Instruction Unit、一個Constant Memory、一個Texture Memory8192Register、一個16KBShare Memory8Stream Processor(SP)和兩個Special Function UnitsSFU)。(GeForce9300M GS只擁有1SM ThreadCUDA模型中最基本的執行單元,執行最基本的程式指令。Block是一組協作ThreadBlock內部允許共享儲存,每個Block最多包含512ThreadGrid是一組Block,共享全域性記憶體。Kernel是在GPU上執行的核心程式,每一個Grid對應一個Kernel任務。 在程式執行的時候,實際上每32Thread組成一個Warp,每個 warp 塊都包含連續的執行緒,遞增執行緒 ID WarpMP的基本排程單位,每次執行的時候,由於MP數量不同,所以一個Block內的所有Thread不一定全部同時執行,但是每個Warp內的所有Thread一定同時執行。因此,我們在定義Block Size的時候應使其為Warp Size的整數倍,也就是Block Size應為32的整數倍。理論上Thread越多,就越能彌補單個Thread讀取資料的latency ,但是當Thread越多,每個Thread可用的暫存器也就越少,嚴重的時候甚至能造成Kernel無法啟動。因此每個Block最少應包含64Thread,一般選擇128或者256,具體視MP數目而定。一個MP最多可以同時執行768Thread,但每個MP最多包含8Block,因此要保持100%利用率,Block數目與其Size有如下幾種設定方式: Ø 2 blocks x 384 threads Ø 3 blocks x 256 threads Ø 4 blocks x 192 threads Ø 6 blocks x 128 threads Ø 8 blocks x 96 threads 

這些電很重要啊,必須要充!不然,我就很難理解為什麼網路執行緒如何分配的。

6 規約思想和同步概念

擴大點說,平行計算是有一種基本思想的,這個演算法能解決很多很常規的問題,而且很實用,比如說累加和累積等——規約思想。對於基礎的、重要的,我想有必要系統的學習。

我覺得有必要重新複製下之前寫的這篇介紹:

http://www.cnblogs.com/viviman/archive/2012/11/21/2780286.html

並行程式的開發有其不同於單核程式的特殊性,演算法是重中之重。根據不同業務設計出不同的並行演算法,直接影響到程式的效率。因此,如何設計並行程式的演算法,似乎成為並程式設計的最大難點。觀其演算法,包括cuda sdk的例子和網上的牛人,給出的一些例子,以矩陣和向量處理為主,深入點的包括fftjulia等數學公式,再高階一點的算是圖形處理方面的例子。學習這些演算法的思想,免不了有自己的一點點總結。之前學習過omp程式設計,結合現在的cuda,我覺得要理解並行程式設計,首先理解劃分和規約這兩個概念。也許你的演算法學的更加紮實。劃分是《演算法》裡面的一個重要思想,將一個大的問題或任務,分解成小問題小任務,各個擊破,最後歸併結果;規約是《cuda**》書上介紹的一個入門的重要思想,規約演算法(reduction)用來求連加、連乘、最值等,應用廣泛。每次迴圈參加運算的執行緒減少一半。不管演算法的思想如何花樣,萬變不離其中的一點--將一個大的任務分解成小的任務集合,分解原則是粒度合適儘量小、資料相關性儘量小。如此而已。因為,我們用GPU是為了加速,要加速必須提高執行任務的並行度!明白這個道理,那麼我們將絞盡腦汁地去想方設法分析自己手上的任務,分解、分解、分解!這裡拿規約來說事情,因為,規約這個東西,似乎可以拿來單做9*9乘法表來熟悉,熟悉了基礎的口訣,那麼99*99的難題也會迎刃而解。ex:向量加法,要實現N=64*256長度的向量的累加和。假設a+b計算一次耗時t

cpu計算:顯然單核的話需要64*256*t。我們容忍不了。

gpu計算:最初的設想,我們如果有個gpu能同時跑N/2個執行緒,我們這N/2個執行緒同時跑,那麼不是隻需要t時間就能將N個數相加程式設計N/2個數相加了嗎?對的。這一輪我們用了t時間;接著的想法,我們不斷的遞迴這個過程,能發現嗎?第二輪,我們用N/2/2個執行緒同時跑,剩下N/2/2個數相加,這一輪我們同樣用了t時間;一直這樣想下去吧,最後一輪,我們用了1個執行緒跑,剩下1個數啦,這就是我們的結果!每一輪時間都為t,那麼理想情況,我們用了多少輪這樣的計算呢?計算次數=log(N)=6*8=48,對的,只用了48輪,也就是說,我們花了48*t的時間!

規約就是這樣,很簡單,很好用,我們且不管程式後期的優化,單從這個演算法分析上來說,從時間複雜度N降到了logN,這在常規演算法上,要提高成這樣的效率,是不得了的,這是指數級別的效率提高!所以,你會發現,GPUCPU無法取代的得天獨厚的優勢——處理單元真心多啊!

規約求和的核函式程式碼如下:

__global__ void RowSum(float* A, float* B)int bid = blockIdx.x; int tid = threadIdx.x; __shared__ s_data[128]; //read data to shared memory  s_data[tid] = A[bid*128 + tid];  __synctheads(); //sync for(int i=64; i>0; i/=2) {  if(tid<i) s_data[tid] = s_data[tid] + s_data[tid+i] ;  __synctheads(); } if(tid==0)   B[bid] = s_data[0];}

這個例子還讓我學到另一個東西——同步!我先不說同步是什麼,你聽我說個故事:我們調遣了10個小組從南京去日本打仗,我們的約定是,10個組可以自己行動,所有組在第三天在上海機場會合,然後一起去日本。這件事情肯定是需要處理的,不能第1組到了上海就先去日本了,這些先到的組,唯一可以做的事情是——等待!這個先來後到的事情,需要統一管理的時候,必須同步一下,在上海這個地方,大家統一下步調,快的組等等慢的組,然後一起幹接下去的旅程。

是不是很好理解,這就是同步在生活中的例子,應該這樣說,計算機的所有機制和演算法很多都是源於生活!結合起來,理解起來會簡單一點。

CUDA中,我們的同步機制用處大嗎?又是如何用的呢?我告訴你,一個正常規模的工程中,一般來說資料都會有先來後到的關係,這一個計算結果可能是提供給另一個執行緒用的,這種依賴關係存在,會造成同步的應用。

__synctheads()這句話的作用是,這個block中的所有執行緒都執行到此的時候,都聽下來,等所有都執行到這個地方的時候,再往下執行。

7 撬開程式設計的鎖

鎖是資料相關性裡面肯定要用到的東西,很好,生活中也一樣,沒鎖,家裡不安全;GPU中沒鎖,資料會被“盜”。

對於存在競爭的資料,CUDA提供了原子操作函式——ATOM操作。

先亮出使用的例子:

__global__ void kernelfun(){__shared__ int i=0;atomicAdd(&i, 1);}

如果沒有加互斥機制,則同一個half warp內的執行緒將對i的操作混淆林亂。

用原子操作函式,可以很簡單的編寫自己的鎖,SKD中有給出的鎖結構體如下:

#ifndef __LOCK_H__#define __LOCK_H__#include "cuda_runtime.h"#include "device_launch_parameters.h"#include "atomic_functions.h"struct Lock {    int *mutex;    Lock( void ) {        HANDLE_ERROR( cudaMalloc( (void**)&mutex, sizeof(int) ) );        HANDLE_ERROR( cudaMemset( mutex, 0, sizeof(int) ) );    }    ~Lock( void ) {        cudaFree( mutex );    }    __device__ void lock( void ) {        while( atomicCAS( mutex, 0, 1 ) != 0 );    }    __device__ void unlock( void ) {        atomicExch( mutex, 0 );    }};#endif


8 CUDA軟體體系結構

 

9 利用好現有的資源

如果連開方運算都需要自己去編寫程式實現,那麼我相信程式設計師這個職業將會縮水,沒有人願意去幹這種活。我想,程式設計師需要學會“偷懶”,現有的資源必須學會高效率的使用。當c++出現了STL庫,c++程式設計師的開發效率可以說倍增,而且程式穩定性更高。

CUDA有提供給我們什麼了嗎?給了,其實給了很多。

先介紹幾個庫:CUFFTCUBLASCUDPP

這裡我先不詳細學習這些庫裡到底有哪些函式,但是,大方向是需要了解的,不然找都不知道去哪兒找。CUFFT是傅立葉變換的庫,CUBLAS提供了基本的矩陣和向量運算,CUDPP提供了常用的並行排序、搜尋等。

CUDA4.0以上,提供了一個類似STL的模板庫,初步窺探,只是一個類似vector的模板型別。有map嗎?map其實是一個散列表,可以用hashtable去實現這項機制。

SDK裡面有很多例子,包括一些通用的基本操作,比如InitCUDA等,都可以固化成函式元件,供新程式的呼叫。

具體的一些可以固化的東西,我將在以後的學習中歸納總結,豐富自己的CUDA庫!

http://blog.csdn.net/huangfengxiao/article/details/8732789

http://blog.csdn.net/huangfengxiao/article/details/8732790

http://blog.csdn.net/huangfengxiao/article/details/8732791

           

給我老師的人工智慧教程打call!http://blog.csdn.net/jiangjunshow

這裡寫圖片描述

5 GPU也不允許偏心

並行的事情多了,我們作為GPU的指令分配者,不能偏心了——給甲做的事情多,而乙沒事做,個麼甲肯定不爽的來。所以,在GPU中,叫做執行緒網路的分配。首先還是來看下GPU的執行緒網路吧,圖2

執行緒網路

我們將具體點的,在主機函式中如果我們分配的是這樣的一個東西:

dim3 blocks(32, 32);

dim3 threads(16, 16);

dim3是神馬?dim3是一個內建的結構體,和linux下定義的執行緒結構體是個類似的意義的東西,dim3結構變數有xyz,表示3維的維度。不理解沒關係,慢慢看。

kernelfun<<<blocks, threads>>>();

我們呼叫kernelfun這個核心函式,將blocksthreads傳到<<<,>>>裡去,這句話可牛逼大了——相當於發號施令,命令那些執行緒去幹活。這裡使用了32*32 * 16*16個執行緒來幹活。你看明白了嗎?blocks表示用了二維的32*32block組,而每個block中又用了16*16的二維的thread組。好吧,我們這個施令動用了262144個執行緒!我們先不管GPU內部是如何排程這些執行緒的,反正我們這一句話就是用了這麼多執行緒。

那我們的核心函式kernelfun()如何知道自己執行的是哪個執行緒?這就是執行緒網路的特點啦,為什麼叫網路,是有講究的,網路就可以定格到網點:

比如int tid = threadId.x + blockId.x * 16

這裡有一個講究,block是有維度的,一維、二維、三維。

對於一維的blocktid = threadId.x

對於(DxDy)二維的blocktid = threadId.x + Dx*threadId.y

對於(DxDyDz)三維的blocktid = threadId.x + Dx*threadId.y + Dz*Dy*threadId.z

我習慣的用這樣的模式去分配,比較通用:

dim3 dimGrid();

dim3 dimBlock();

kerneladd<<<dimGrid, dimBlock>>>();

這可是萬金油啊,你需要做的事情是填充dimGriddimBlock的結構體建構函式變數,比如,dimGrid(16, 16)表示用了16*16的二維的block執行緒塊。

0,0)(0,1)(0,2)……(0,15

1,0)(1,1)(1,2)……(1,15

2,0)(2,1)(2,2)……