来源：http://luofl1992.is-programmer.com/posts/38847.html
编程语言的特点是要实践，实践多了才有经验。很多东西书本上讲得不慎清楚，不妨自己用代码实现一下。
作为例子，我参考了书本上的矩阵相乘的例子，这样开始写代码，然后很自然地出现了各种问题。
以下的内容供大家学习参考，有问题可以留言与我反馈。
开始学着使用 CUDA，实现一个矩阵乘法运算。
首先我们要定义一个矩阵的结构体，话说CUDA是否支持结构体作为设备端的函数的参数呢？
不妨都一股脑试验一下。
1、安装CUDA 5.0
在NVIDIA CUDA官方网站下载对应自己操作系统的最新版 CUDA，
为什么要最新版本呢，一方面技术一直更新，自己要保持潮流嘛~（程序员的苦逼之处也在此）
更新的版本一般实现功能上更加完善，BUG更少。
废话不多说，去这里下载：
https://developer.nvidia.com/cuda-downloads
最新版的CUDA不像4.0版本分ToolKit、SDK、Samples，可以说简化了安装复杂度。
安装完毕可以看看示例，路径为
安装目录 \\ nvToolsExt\samples
比如我的为：C:\Program Files\NVIDIA GPU Computing Toolkit\nvToolsExt\samples
2、创建一个新的CUDA项目
CUDA 5.0不需要再额外安装 cudaVSWizard了，安装完成后，VS中会出现这样的东西，
新建 ===》 项目 ===》 NVIDIA =====》 CUDA
填写项目名称，然后就会自动生成一个kernal.cu文件，编辑这个文件就可以了。
如果代码比较复杂，可以弄多个文件，但是cuda目前所有的设备代码似乎必须写在一个源文件中，不能使用常用的函数声明+其他文件实现的方式。如果存在多个cu文件，可以用#include "axxx.cu"指令打包成一个cu源文件。
这个后续还要测试。贴一张图说明一下，我的为VS2010：
3、修改代码
关于矩阵运算，比例子复杂一些。这里先给出其核心代码，
相信有良好C语言基础的人能够轻松看懂这个函数而不需要我的注释。
?
矩阵乘法的CUDA C代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
|
// 矩阵乘法，C = A * B
__global__ void gpuMatrixMul(const Matrix A, const Matrix B, Matrix C)
{
int row = blockIdx.x * blockDim.x + threadIdx.x;
int col = blockIdx.y * blockDim.y + threadIdx.y;
float ret = 0;
if ( row >= A.row || col >= B.col || A.col != B.row )
return;
float *ma = A.ele + row * A.stride;
float *mb = B.ele + col;
for ( int i = 0, n = A.col, stride = B.stride; i < n; i++ )
{
ret += ma[ i ] * mb[ i * stride ];
}
C.ele[row * C.stride + col] = ret;
}
/*
__global__ void gpuMatrixMul(const Matrix *A, const Matrix *B, Matrix *C)
{
int row = blockIdx.x * blockDim.x + threadIdx.x;
int col = blockIdx.y * blockDim.y + threadIdx.y;
float ret = 0;
if ( row >= A->row || col >= B->col || A->col != B->row )
return;
float *ma = A->ele + row * A->stride;
float *mb = B->ele + col;
for ( int i = 0, n = A->col, stride = B->stride; i < n; i++ )
{
ret += ma[ i ] * mb[ i * stride ];
}
C->ele[row * C->stride + col] = ret;
}
*/
|
其中结构体Matrix的定义如下：
?
Matrix结构体定义
1
2
3
4
5
6
7
8
|
// 矩阵结构体
typedef struct MartixTag
{
float *ele; // 一维数组，保存矩阵的元素
size_t col; // 矩阵的列数
size_t row; // 矩阵的行数
size_t stride; // 矩阵一行数据的数量，为了方便内存对齐，2^n
}Matrix;
|
实际测试表明，上述代码中，采用指针传递结构体参数时，运算结果会不正确。
这点在我们以后写CUDA函数的时候应该注意，我的理解是这样子，不知道大家怎么看：
?
CUDA函数参数注意
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
|
// （1）设备端分配的显存空间，不能直接在CPU上访问；
// （2）主机端分配的内存空间，不能在设备端使用；
// （3）需要使用时，需要使用 CUDA的相关函数。
cudaError_t cudaMalloc(void **p, size_t size);
/*
这个cudaMalloc函数中 参数p指针的指针，
主要是为了能够修改改指针指向的地址，而函数返回类型可以指出是否出现了错误。
这个部分可以再回顾一下指针的概念。
*/
cudaError_t cudaMemcpy(void *dst, void *src, size_t cnt, cudaMemoryKind kind);
/*
这个cudaMemcpy函数中
cnt指定要复制的内存字节数。
kind指定了拷贝的方向，主要有
cudaMemcpyHostToDevice // CPU内存到GPU显存
cudaMemcpyDeviceToHost // GPU显存到CPU内存
*/
|
4、调试设备端代码
在VS平台上使用，建议安装Nsight，这在官方网站上下载就好了。
需要自己注册一个NVIDIA帐号，然后填写一些信息，
什么？页面全英文看不懂？
慢慢看吧。https://developer.nvidia.com/nvidia-nsight-visual-studio-edition
显卡驱动，一般以前安装很新的版本的，就不用再安装了。
安装一下Nsight，关闭VS，完毕后重新打开时，发现菜单栏里面多了一个选项。
在要调试的设备端代码中按F9设置一个断点，然后选中菜单栏：
Nsight ---> Start CUDA Debuging ---->
就可以看到代码会中断在设备端的代码中了，设置在CPU端代码的断点不会被识别。
5、附录代码
最后附上我的一个成果，按照书上的分块计算矩阵乘法代码，经本人测试无误：
?
分块计算矩阵乘积
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
|
__global__ void gpuMatrixMulB(Matrix A, Matrix B, Matrix C)
{
if ( blockDim.x != BLOCK_SIZE || blockDim.y != BLOCK_SIZE )
return;
// x 横向指向一行 // y 指向列号，数组的第二维
// 访问元素以 C.ele[x][y] 的形式
// 一个thread计算C的一个值
const int bx = blockIdx.x, by = blockIdx.y;
const int tx = threadIdx.x, ty = threadIdx.y;
// A中第一个子块的起始地址
const int aBegin = A.stride * bx * BLOCK_SIZE;
const int aEnd = aBegin + A.col;
const int aStep = BLOCK_SIZE;
// B中要处理的第一个子块的起始地址
const int bBegin = BLOCK_SIZE * by;
const int bEnd = bBegin + B.stride * B.row;
const int bStep = BLOCK_SIZE * B.stride;
float Csub = 0;
// 循环A的一行（by指向最开始），和B的一列（bx指向最开始）
for ( int a = aBegin, b = bBegin, i = 0; a < aEnd; a += aStep, b+= bStep, i++ )
{
__shared__ float AS[BLOCK_SIZE][BLOCK_SIZE];
__shared__ float BS[BLOCK_SIZE][BLOCK_SIZE];
int na = min(BLOCK_SIZE, aEnd - a);
// 对不足一个整块的不用整体复制
if ( tx + bx * BLOCK_SIZE < A.row && ty + a - aBegin < A.col )
AS[tx][ty] = A.ele[a + A.stride * tx + ty];
else
AS[ty][tx] = 0;
if ( tx * B.stride + b < bEnd && ty + by * BLOCK_SIZE < B.col )
BS[tx][ty] = B.ele[b + B.stride * tx + ty];
else
BS[tx][ty] = 0;
__syncthreads();
for ( int k = 0 ; k < na; ++k )
{
Csub += AS[tx][k] * BS[k][ty];
}
__syncthreads();
}
if ( by * BLOCK_SIZE + ty > C.col || bx * BLOCK_SIZE + tx > C.row )
return;
// C[by*BLOCK + ty][bx * BLOCK + tx]
C.ele[ (by * BLOCK_SIZE + ty ) + C.stride * (bx*BLOCK_SIZE + tx) ] = Csub;
}
|
调用方式如下（对相关参数初始化完毕后执行）：
?
调用方式
1
2
3
4
5
6
7
8
9
10
11
12
13
14
|
// 启用 kernel
dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
dim3 dimGrid( (d_C.row + dimBlock.x - 1) / dimBlock.x, (d_C.col + dimBlock.y - 1) / dimBlock.y);
fprintf( stderr, "DimBlock: %d, %d, %d\nDimGrid: %d, %d, %d\n",
dimBlock.x, dimBlock.y, dimBlock.z, dimGrid.x, dimGrid.y, dimGrid.z);
#if 1
// 不要通过指针参数调用，指针所指向的内存必须要通过显式的内存拷贝
// 才能在GPU中使用
// gpuMatrixMul<<<dimGrid, dimBlock>>>(&d_A, &d_B, &d_C);
gpuMatrixMulB<<<dimGrid, dimBlock>>>( d_A, d_B, d_C);
#else
gpuMatrixMul<<<dimGrid, dimBlock>>>( d_A, d_B, d_C);
#endif
cudaMemcpy( c.ele, d_C.ele, c.stride * c.row * sizeof(float), cudaMemcpyDeviceToHost );
|
顺便附录一份C源代码计算矩阵乘积，也分享给大家：
?
C语言计算矩阵乘积
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
|
int MatMul(const Matrix *a, const Matrix *b, Matrix *pRet)
{
const size_t row = a->row;
const size_t col = b->col;
const size_t n = a->col;
const size_t stride = b->stride;
float const * ma = a->ele;
float const * mb = b->ele;
float * mRet = pRet->ele;
int i = 0, j = 0, k = 0;
float ret = 0;
if ( a->col != b->row )
return 0;
if ( pRet->ele == NULL )
{
pRet->col = a->col;
pRet->row = b->row;
pRet->stride = a->stride;
pRet->ele = (float *)malloc( a->stride * sizeof(float) * b->row );
}
for ( i = 0; i < row; i++ )
{
for ( j = 0; j < col; j++ )
{
ret = 0.0;
for ( k = 0; k < n; k++ )
{
ret += ma[k] * mb[k * stride + j ];
}
mRet[j] = ret;
}
ma += a->stride;
mRet += pRet->stride;
}
return 1;
}
|
时间： 2024-10-06 06:39:38
