(2)那么下面就是不使用shared
memory的并行化算法的思路。简单地来说,就是将上述可并行化的部分传递给GPU,使用CUDA来计算。代码如下:
void MatrixMulOnDevice(float*
M, float* N, float* P, intWidth)
{
int size = Width * Width * sizeof(float);
float* Md,
Nd,
Pd;
//设置调用内核函数时的线程数目
dim3 dimBlock(Width,
Width);
dim3 dimGrid(1,
1);
//在设备存储器上给M和N矩阵分配空间,并将数据复制到设备存储器中
cudaMalloc(&Md,
size);
cudaMemcpy(Md,
M, size, cudaMemcpyHostToDevice);
cudaMalloc(&Nd,
size);
cudaMemcpy(Nd,
N, size, cudaMemcpyHostToDevice);
//在设备存储器上给P矩阵分配空间
cudaMalloc(&Pd,
size);
//内核函数调用,将在后续部分说明
//只使用了一个线程块(dimGrid),此线程块中有Width*Width个线程
MatrixMulKernel<<
dimBlock>>>(Md,
Nd,
Pd,Width);
// 从设备中读取P矩阵的数据
cudaMemcpy(P,
Pd, size, cudaMemcpyDeviceToHost);
// 释放设备存储器中的空间
cudaFree(Md);
cudaFree(Nd);
cudaFree (Pd);
}
__global__ void MatrixMulKernel(float*
Md,
float* Nd,
float* Pd, intWidth)
{
// 2维的线程ID号
int tx = threadIdx.x;
int ty = threadIdx.y;
// Pvalue用来保存被每个线程计算完成后的矩阵的元素
float Pvalue = 0;
//每个线程计算一个元素
for
(intk
= 0; k < Width; ++k)
{
float
Melement = Md[ty
* Width + k];
float
Nelement = Nd[k
* Width + tx];
Pvalue
+= Melement * Nelement;
}
// 将计算结果写入设备存储器中
Pd[ty*
Width + tx]
= Pvalue;
}
(3)这里面有一个很重要的问题就是所有的数据都放在了global
memory里面,数据的传输速度会受到带宽的限制,并且矩阵的大小受到计算能力的限制。所以可以将存储器优化,算法优化,使它更有效率地执行。
思路如下:
将整个大的width*width的矩阵分成tile_width*tile_width的小矩阵,一共可分为(width/tile_width)*(width/tile_width)个。当然此处要注意tile_width和计算能力的大小。
然后是将这些小矩阵分别进行计算,将计算结果放入shared
memory中,如上面的第二张图。简单来说就是M的蓝色部分乘以N的蓝色部分,结果放入黄色部分中;再用M的橙色部分乘以N的橙色部分,结果加入黄色部分中,依次类推。
定义blocksPerGrid和threadsPerBlock:
dim3
dimBlock(tile_width,tile_width);
dim3
dimGrid(width/tile_width,width/tile_width);
调用内核函数:
MatrixMulKernel<<>>(Md,Nd,Pd,width);
内核函数如下:
__global__ void MatrixMulKernel(float*
Md,
float* Nd,
float* Pd, intWidth)
{
//获得线程块号
intbx=
blockIdx.x;
intby
= blockIdx.y;
//获得块内的线程号
inttx=
threadIdx.x;
intty=
threadIdx.y;
//Pvalue:线程计算完成后的子矩阵元素——自动变量
float Pvalue = 0;
//循环,遍历M和N的所有子矩阵
for (intm
= 0; m < Width/TILE_WIDTH; ++m)
{
//获取指向当前矩阵M子矩阵的指针Msub
Float* Mdsub =
GetSubMatrix(Md,m,by,Width);
//获取指向当前矩阵N的子矩阵的指针Nsub
Float* Ndsub = GetSubMatrix(Nd,bx,m,Width);
//共享存储器空间声明
__shared__ float Mds[tile_width][tile_width];
__shared__ float Nds{tile_width][tile_width];
//每个线程载入M的子矩阵的一个元素
Mds[ty][tx] =
GetMatrixElement(Mdsub,tx,ty);
//每个线程载入N的子矩阵的一个元素
Nds[ty][tx] = GetMatrixElement(Ndsub,tx,ty);
//同步,在计算之前,确保子矩阵所有的元素都已经载入共享存储器中
__syncthreads();
//每个线程计算线程块内子矩阵中的一个元素
for(int
k=0;k
Pvalue+=Mds[ty][k]*Nds[k][tx];
//同步,确保重新载入新的M和N子矩阵数据前,上述计算过程已经全部完成
__syncthreads();
};
}
以下是GetSubMatrix(Md,x,y,Width)【获取第(x,y)号子矩阵的起始地址】{Md+y*tile_width*width+x*tile_width;}和GetMatrixElement(Mdsub,tx,ty,Width)【获取子矩阵中某个元素的地址】{*(Mdsub+ty*width+tx);}的示意图:
接下来是保存结果:
//获取指向矩阵P的子矩阵的指针
Matrix Psub
= GetSubMatrix(P,bx,by);
//向全局存储器写入线程块计算后的结果子矩阵
//每个线程写入一个元素
SetMatrixElement(Psub,tx,ty,Pvalue);
总的来说,内核函数的思路就是:
①获得线程块号以及块内的线程号
②申明一个自动变量以记录线程计算完成之后的子矩阵元素
③从第一个子矩阵开始,遍历M和N的所有子矩阵,以完成如下操作:
④获取指向当前矩阵M、N子矩阵的指针Msub和Nsub
⑤声明两个tile_width*tile_width大小的共享存储器空间,分别对应M的子矩阵和N的子矩阵
⑥分别载入M和N的子矩阵的元素到共享存储器中
⑦同步,确保所有的元素都已经载入共享存储器中
⑧开始计算,每个线程计算线程块内子矩阵的一个元素
⑨同步,确保计算过程已经完成
⑩获取指向矩阵P的子矩阵的指针
⑪向全局存储器中写入计算结果后的子矩阵
但是小白现在还有一个问题,按理说,Pvalue应该是存储在register中才对,那么,向全局存储器中写入计算结果的值是CUDA自动从寄存器中取值然后再写入吗,但又有听说register是不可寻址的?所以表示不清楚,只好留待以后探讨吧。