這幾天研究了一下CUDA，發(fā)現(xiàn)其并行的思想和普通的CPU多線程思想不太一致，但還是挺不錯。主要是將任務(wù)劃分成一個個block，然后每個block里面再劃分成細的線程。然后每個線程做自己做的
事情。這種并行思想很適用于像矩陣運算這些元素與元素之間的運算并不耦合得很厲害，但整體數(shù)據(jù)很大的情況，這只是我對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;耦合程度很小，所以我們可以通過劃分區(qū)域的方法，讓每個線程負責一個區(qū)域。
怎么劃分呢？首先最初的想法是讓每一個線程計算一個C(i,j)，那么估算一下，應(yīng)該需要height_c*width_c，也就是ha*wb個線程。進一步，我們將矩陣按一個大方格Grid劃分，如果一個
方格Grid大小是16*16，那么矩陣80*48的可以表示為5(*16) * 3(*16)，即16*16個大格子(block)，每一個格子內(nèi)，自然就是(height_c/16) *(width_c/16)個線程了。
好了，劃分完后，內(nèi)核代碼如下：
計算版本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) }發(fā)現(xiàn)，我們還可以再細度劃分：
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的子乘積，計算完畢后加起來則可。這里主要使用了共享顯存作優(yōu)化。
計算版本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的接口函數(shù)：
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);
}
調(diào)用的主函數(shù)如下：
#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;
}
運算結(jié)果如下：
版本0：
版本1：
可以看出，GPU并行性能比CPU好很多，而且版本1優(yōu)于版本0
整個工程下載：
/Files/bennycen/CUDAMatrixMul.rar
posted on 2011-07-26 17:01
bennycen 閱讀(4641)
評論(1) 編輯 收藏 引用 所屬分類:
CUDA