(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是不可寻址的？所以表示不清楚，只好留待以后探讨吧。