下面的例子代码利用了共享存储器实现矩阵相乘。本实现中，每个线程块负责计算一个小方阵Csub，Csub是C的一部分，而块内的每个线程计算Csub的一个元素。如图3-2所示。Csub等于两个长方形矩阵的乘积：A的子矩阵尺寸是(A.width,block_size)，行索引与Csub相同，B的子矩阵的尺寸是(block_size,A.width)，列索引与Csub相同。为了满足设备的资源，两个长方形的子矩阵分割为尺寸为block_size的方阵，Csub是这些方阵积的和。每次乘法的计算是这样的，首先从全局存储器中将二个对应的方阵载入共享存储器中，载入的方式是一个线程载入一个矩阵元素，然后一个线程计算乘积的一个元素。每个线程积累每次乘法的结果并写入寄存器中，结束后，再写入全局存储器。
采用这种将计算分块的方式，利用了快速的共享存储器，节约了许多全局存储器带宽，因为在全局存储器中，A只被读了(B.width/block_size)次同时B读了(A.height/block_size)次。
前面代码中的Matrix 类型增加了一个stride域，这样子矩阵能够用同样的类型有效表示。__device__函数(见B.1.1节)用于读写元素和从矩阵中建立子矩阵。
// Matrices are stored in row-major order:
// M(row, col) = *(M.elements + row * M.stride + col)
typedef struct {
int width;
int height;
int stride;
float* elements;
} Matrix;
// Get a matrix element
__device__ float GetElement(const Matrix A, int row, int col) {
return
// M(row, col) = *(M.elements + row * M.stride + col)
typedef struct {
int width;
int height;
int stride;
float* elements;
} Matrix;
// Get a matrix element
__device__ float GetElement(const Matrix A, int row, int col) {
return
// Read C from device memory
cudaMemcpy(C.elements, d_C.elements, size,
cudaMemcpyDeviceToHost);
// Free device memory
cudaFree(d_A.elements);
cudaFree(d_B.elements);
cudaFree(d_C.elements);
}
cudaMemcpy(C.elements, d_C.elements, size,
cudaMemcpyDeviceToHost);
// Free device memory
cudaFree(d_A.elements);
cudaFree(d_B.elements);
cudaFree(d_C.elements);
}
