這幾天研究了一下CUDA，發現其並行的思想和普通的CPU多線程思想不太一致，但還是挺不錯。主要是將任務劃分成一個個block，然後每個block裡面再劃分成細的線程。然後每個線程做自己做的
事情。這種並行思想很適用於像矩陣運算這些元素與元素之間的運算並不耦合得很厲害，但整體數據很大的情況，這只是我對CUDA的初步感覺。
矩陣相乘的CPU程序如下：
//C = A*B
void MatrixMulCPU(float* _C,const float *_A,const float *_B,int _wa,int _ha,int _wb)
{
float sum = 0;
for (int i = 0; i < _ha; ++i)
{
for (int j = 0; j < _wb; ++j)
{
sum = 0;
for (int k = 0; k < _wa; ++k)
{
sum += (float)_A[i*_wa+k]*(float)_B[k*_wb+ j];
}
_C[i*_wb+j] = (float)sum;
}
}
}
從上面可以看出，C(i,j) = sum { A(i,k)*B(k,j) } 0<=k < _wa;耦合程度很小，所以我們可以通過劃分區域的方法，讓每個線程負責一個區域。
怎麼劃分呢？首先最初的想法是讓每一個線程計算一個C(i,j)，那麼估算一下，應該需要height_c*width_c，也就是ha*wb個線程。進一步，我們將矩陣按一個大方格Grid劃分，如果一個
方格Grid大小是16*16，那麼矩陣80*48的可以表示為5(*16) * 3(*16)，即16*16個大格子(block)，每一個格子內，自然就是(height_c/16) *(width_c/16)個線程了。
好了，劃分完後，內核代碼如下：
計算版本0：
__global__ void matrix_kernel_0(float* _C,const float* _A,const float *_B,int _wa,int _wb)
{
float sum = 0;
//找出該線程所在的行列
int row = blockIdx.y*blockDim.y + threadIdx.y;
int col = blockIdx.x*blockDim.x + threadIdx.x;
//線程Thread(row,col)負責計算C(row,col)
for (int i = 0; i < _wa; ++i)
{
sum += _A[row*_wa + i]*_B[i*_wb + col];
}
_C[row*_wb + col] = sum;
}
另外一種思路，我們不讓每一個線程完整計算一個C(i,j)，通過C(i,j) = sum { A(i,k)*B(k,j) }發現，我們還可以再細度劃分：
Csub(i,j) = sum{A(i,ksub+offsetA)*B(ksub+offsetB,j)} 0<=ksub < blockSize
C(i,j) = sum{Csub(i,j)}
就是把矩陣分成n*n個大的子塊，然後每一個block負責計算子塊i 和 子塊j的子乘積，計算完畢後加起來則可。這裡主要使用了共享顯存作優化。
計算版本1：
__global__ void matrix_kernel_1(float* _C,const float* _A,const float *_B,int _wa,int _wb)
{
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
//該block要處理的A
int aBegin = _wa*(by*BLOCK_SIZE);//A(0,by)
int aEnd = aBegin + _wa - 1;
int aStep = BLOCK_SIZE;//offsetA
int bBegin = BLOCK_SIZE*bx;//B(bx,0)
int bStep = BLOCK_SIZE*_wb;//offsetB
float cSub = 0;
for (int a = aBegin,b = bBegin; a <= aEnd; a += aStep,b += bStep)
{
__shared__ float As[BLOCK_SIZE][BLOCK_SIZE];
__shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];
//每個線程負責一個元素拷貝
As[ty][tx] = _A[a + _wa*ty + tx];
Bs[ty][tx] = _B[b + _wb*ty + tx];
__syncthreads();
//每個線程負責計算一個子塊i 和 子塊j的子乘積
for (int k = 0; k < BLOCK_SIZE; ++k)
{
cSub += As[ty][k]*Bs[k][tx];
}
__syncthreads();
}
//全局地址，向全局寄存器寫回去
//一個線程負責一個元素，一個block負責一個子塊
int cIndex = (by*BLOCK_SIZE + ty)*_wb + (bx*BLOCK_SIZE + tx);
_C[cIndex] = cSub;
}
最後寫一個面向Host的接口函數：
void matrixMulGPU(float* _C,const float *_A,const float *_B,int _wa,int _ha,int _wb)
{
float* d_a = myNewOnGPU<float>(_wa*_ha);
float* d_b = myNewOnGPU<float>(_wb*_wa);
float* d_c = myNewOnGPU<float>(_wb*_ha);
copyFromCPUToGPU(_A,d_a,_wa*_ha);
copyFromCPUToGPU(_B,d_b,_wb*_wa);
dim3 threads(BLOCK_SIZE,BLOCK_SIZE);
dim3 blocks(WC/BLOCK_SIZE,HC/BLOCK_SIZE);
matrix_kernel_0<<<blocks,threads>>>(d_c,d_a,d_b,_wa,_wb);
cudaThreadSynchronize();
copyFromGPUToCPU(d_c,_C,_wb*_ha);
myDeleteOnGPU(d_a);
myDeleteOnGPU(d_b);
myDeleteOnGPU(d_c);
}
調用的主函數如下：
#include <stdio.h>
#include <cuda_runtime.h>
#include <cutil.h>
#include <cutil_inline.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
#include <Windows.h>
#include "CUDACommon.h"
#include "MatrixMulCPU.h"
#include "MatrixMulGPU.h"
void randomInit(float* _data,int _size)
{
for (int i = 0; i < _size; ++i)
{
_data[i] = rand()/(float)RAND_MAX;
}
}
bool checkError(const float* _A,const float* _B,int _size)
{
for (int i = 0 ; i < _size; ++i)
{
if (fabs(_A[i] - _B[i]) > 1.0e-3)
{
printf("%f \t %f\n",_A[i],_B[i]);
return false;
}
}
return true;
}
int main(int argc, char* argv[])
{
srand(13);
if(!InitCUDA()) {
return 0;
}
float* A = myNewOnCPU<float>(WA*HA);
float* B = myNewOnCPU<float>(WB*HB);
randomInit(A,WA*HA);
randomInit(B,WB*HB);
float* C = myNewOnCPU<float>(WC*HC);
memset(C,0,sizeof(float)*WC*HC);
float* C2 = myNewOnCPU<float>(WC*HC);
memset(C2,0,sizeof(float)*WC*HC);
unsigned int tick1 = GetTickCount();
MatrixMulCPU(C2,A,B,WA,HA,WB);
printf("CPU use Time : %dms\n",GetTickCount() - tick1);
unsigned int timer = 0;
cutilCheckError(cutCreateTimer(&timer));
cutilCheckError(cutStartTimer(timer));
{
matrixMulGPU(C,A,B,WA,HA,WB);
}
cutilCheckError(cutStopTimer(timer));
printf("GPU use time: %f (ms) \n", cutGetTimerValue(timer));
cutilCheckError(cutDeleteTimer(timer));
if (checkError(C,C2,WC*HC))
{
printf("Accept\n");
}
else
{
printf("Worng Answer\n");
}
myDeleteOnCPU(A);
myDeleteOnCPU(B);
myDeleteOnCPU(C);
myDeleteOnCPU(C2);
return 0;
}
運算結果如下：
版本0：
版本1：
可以看出，GPU並行性能比CPU好很多，而且版本1優於版本0
