前言
GPU有許多不同的記憶體種類(global memory, shared memory, texture memory, register, local memory, etc..)，其中最常用的是global memory，而global memory的使用在xxx(還沒寫＠＠)也都有舉例說明如何去使用，一般的流程是，先配置好global
memory，然後再將CPU的資料給複製到device 的global memory中，之後在GPU中提取global memory的值來進行運算並把結果存回global memory中，但是對於某些問題來講，記憶體的存取頻繁，也就有可能因為global memory的使用而造成整體運算效能的降低，因為global memory的存取速度是GPU中的記憶體裡最低的，若我們可以善加利用其他較快的記憶體進行存取，則總體效能應該可以再提升，為了要讓大家了解如何將原本使用global memory存取運算的CUDA程式，改寫成使用shared memory存取運算的CUDA程式，進一步獲得相對應的效能提升，這邊就先以矩陣相乘做範例，來帶大家一步步做喔～
基本的矩陣相乘
矩陣相乘一直是許多科學運算中會去使用到的基本運算問題，
假設A為m*l矩陣(A = [aij] m*l)，B為l*n矩陣(B = [bij] l*n)，則他們的乘積AB會為一個m*n的矩陣C (C = [cij] m*n)，其中
Ci,j = 
，我們可以看成C矩陣的每一個元素就是A的row向量與B的column向量進行向量內積，也就是每一個C的element都可以看成獨立的向量內積運算的問題，這樣的問題很適合來進行平行運算，而過去在CPU中的平行運算架構中已經有不少對於矩陣相乘的研究，而現在讓我們來看看如何利用CUDA改寫在CPU上的矩陣相乘到GPU的運算環境上。
Example 1: CPU版本
//-------------------------------------------------------------------- //A : m x l //B : l x n //C : m x n (C=A*B) //-------------------------------------------------------------------- void host_mm(float* C, float* A, float* B, int m, int n, int l) { for(int i=0; i<m; i++) for(int j=0; j<n; j++) { float s=0; for (int k=0; k<l; k++) { float a = A[i*l + k]; float b = B[k*n + j]; s += a * b; } C[i*n + j] = s; } }
SBMT(Single Block, Multiple Threads)
在Example 1裡的host_mm 函式利用C語言在CPU裡實作了傳統的矩陣相乘運算，我們可以看到他的運算需要使用到三個迴圈，也就是需要O(n^3)的時間複雜度，運算起來很耗費時間，在這函式中我們利用i*l + k 與 k*n + j 的index分別把A與B陣列裡的row與column向量取出進行向量內積的運算，之後將其結果存回C陣列i*n + j的index之中，事實上，這樣的運算我們可以將其分為許多向量內積的問題，一開始我們利用CUDA解決這個問題時，先不要想說要利用多少個thread block來進行，考量最簡單的case，一個thread block的情形，那一個thread block裡有許多thread，要如何分配這些thread去做相對應的運算呢? Figure 1提供了一種方法：
Figure 1: SBMT(Single Block, Multiple Threads) method
矩陣相乘得出的每一個element的運算，可以對應到thread block的其中一個thread的運算，也就是這個thread可以將A,B矩陣的row與column給取出來進行向量內積，得出的結果存到此thread對應的C矩陣的某位址中，以圖3.1.1的例子而言，編號為(2,2)的thrad負責C矩陣(2,2)位址的值，而其中的運算就是將A的第三個row向量(1,1,1,1)與B的第三個column向量(1,2,3,4)，進行向量內積的運算，所以為1*1+1*2+1*3+1*4 = 10，這樣的方法我們可以嘗試著寫一個簡單的CUDA程式看看：
Example 2: (SBMT_matrix_global.cu)
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <sys/time.h>
//--------------------------------------------------------------------
//A : m x l
//B : l x n
//C : m x n (C=A*B)
//--------------------------------------------------------------------
void host_mm(float* C, float* A, float* B, int m, int n, int l){
for(int i=0; i<m; i++)
for(int j=0; j<n; j++)
{
float s=0;
for (int k=0; k<l; k++)
{
float a = A[i*l + k];
float b = B[k*n + j];
s += a * b;
}
C[i*n + j] = s;
}
}
//--------------------------------------------------------------------
__global__ void gpu_mm(float* C, float* A, float* B, int m, int n, int l){
//// 2D Thread ID
int tx = threadIdx.x;
int ty = threadIdx.y;
// Pvalue is used to store the element of the matrix
// that is computed by the thread
float Pvalue = 0;
for (int k = 0; k < l; ++k)
{
float Aelement = A[ty * l + k];
float Belement = B[k * n + tx];
Pvalue += Aelement * Belement;
}
// Write the matrix to device memory;
// each thread writes one element
C[ty * n + tx] = Pvalue;
}
//----------------------------------------------
double diff(float* a, float* b, int n){
double s=0, r=0;
for(int k=0; k<n; k++)
{
double w=a[k]-b[k];
s+=w*w;
r+=a[k]*a[k];
}
return sqrt(s/r);
}
void random(float* a, int n){
for(int k=0; k<n; k++){
a[k]=(float)rand()/RAND_MAX*2-1;
}
}
//----------------------------------------------
void testMatrix(int m, int n, int l)
{
//initialize
float *a = (float*)malloc(sizeof(float)*m*l);
float *b = (float*)malloc(sizeof(float)*l*n);
float *c1 = (float*)malloc(sizeof(float)*m*n);
float *c2 = (float*)malloc(sizeof(float)*m*n);
srand(time(0));
random(a,m*l);
random(b,l*n);
memset(c1, 0, sizeof(float)*m*n);
memset(c2, 0, sizeof(float)*m*n);
float *ga,*gb,*gc;
cudaMalloc((void**)&ga, m*l*sizeof(float));
cudaMalloc((void**)&gb, l*n*sizeof(float));
cudaMalloc((void**)&gc, m*n*sizeof(float));
cudaMemcpy(ga, a, m*l*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(gb, b, l*n*sizeof(float), cudaMemcpyHostToDevice);
cudaMemset(gc, 0, m*n*sizeof(float));
cudaEvent_t start,stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start,0);
//SBMT(Single Block, Multiple Threads)
gpu_mm<<<dim3(1,1,1), dim3(m, n, 1)>>> (gc,ga,gb,m,n,l);
cudaThreadSynchronize();
cudaEventRecord(stop,0);
cudaEventSynchronize(stop);
float elapsedTime;
cudaEventElapsedTime(&elapsedTime, start, stop);
cudaMemcpy(c2, gc, m*n*sizeof(float), cudaMemcpyDeviceToHost);
double c_start,c_stop;
double CPU_execution_time;
c_start = (double)clock();
host_mm(c1, a, b, m, n, l);
c_stop = (double)clock();
CPU_execution_time = (c_stop - c_start)/(double)CLOCKS_PER_SEC;
//check precision
double err=diff(c1,c2,m*n);
printf("err = %g\n", err);
printf(" ======== Execution Infomation ========\n");
printf(" Excuetion Time on GPU: %3.20f s\n",elapsedTime/1000);
printf(" Excuetion Time on CPU: %3.20f s\n",CPU_execution_time);
printf(" Speed up = %f\n",(CPU_execution_time/(elapsedTime/1000)));
printf(" ======================================\n\n");
free(a);
free(b);
free(c1);
free(c2);
cudaFree(ga);
cudaFree(gb);
cudaFree(gc);
}
//----------------------------------------------
int main()
{
int m=32;
int n=32;
int l=32;
testMatrix(m,n,l);
return 0;
}
為kernel launch，dim3為CUDA內建資料結構，就代表配置grid, block的維度，
我們舉的例子分別代表擁有m*n個thread的thread block與擁有1*1個thread block的kernel，dim3(1,1,1), dim3(m, n, 1)，放入"<<< >>>" 這tag之中，這樣便代表呼叫此kernel function的grid擁有1*1個thread block，以及每個block擁有m*n個thread。
kernel function使用了dim3這個內建的資料結構來配置了grid和block的大小，於是我們可以在kernel function裡使用多維度的thread block與thread ID，我們的例子thread block數量只有1，所以沒有多維度的使用問題，而在一個thread block裡的thread數量有m*n個，所以m*n 個thread 都有一個唯一的2D編號來表示，以下程式碼由
int tx = threadIdx.x;int ty = threadIdx.y;
來得到2D thread ID。
在kernel function呼叫之前的cudaMalloc已經配置好global memory大小的陣列了，所以我們在gpu_mm這個kernel function裡面可以直接提取A,B陣列在global memory裡的值，而我們可以發現gpu_mm 函式裡面只要一個for迴圈來抓取A,B矩陣的row 與column值，之後算出來的值再存回C矩陣對應的元素裡。
__global__ void gpu_mm(float* C, float* A, float* B, int m, int n, int l){//// 2D Thread IDint tx = threadIdx.x;int ty = threadIdx.y;// Pvalue is used to store the element of the matrix// that is computed by the threadfloat Pvalue = 0;for (int k = 0; k < l; ++k){float Aelement = A[ty * l + k];float Belement = B[k * n + tx];Pvalue += Aelement * Belement;}// Write the matrix to device memory;// each thread writes one elementC[ty * n + tx] = Pvalue;}
結論
我們的矩陣大小為32*32 = 1024，為了剛好符合我的device一個block的thread個數上限，測試的GPU為GeForce GT 650M，CC(CUDA Capability)為3.0，CPU為2.6 GHz Intel Core-i7，往後沒有特別註明的話都是在此環境下測試，測試的結果如下：
err = 9.88651e-08我們看到GPU的表現並沒有很好，寫了一大堆結果跑得比CPU慢！那ㄟ安ㄋ？別緊張，我們來看一下這數據背後的意義，首先這個時間的數據都很小，因為一開始我們只有用32*32這麼大的矩陣相乘，說實在的執行時間都極短，測量時間的誤差也大，再來就是我們總共只有使用1024個thread，且只有一個block，那同時間也只有一個SM(Streaming Multiprocessor)在運行而已，總之就是沒有充分使用到GPU資源囉，我們接下來再來看要怎麼改善 (待續...)
======== Execution Infomation ========
Excuetion Time on GPU: 0.00069228798383846879 s
Excuetion Time on CPU: 0.00012000000000000000 s
Speed up = 0.173338
======================================
