CUDA学习笔记(二)【转】

来源:http://luofl1992.is-programmer.com/posts/38847.html

编程语言的特点是要实践,实践多了才有经验。很多东西书本上讲得不慎清楚,不妨自己用代码实现一下。

作为例子,我参考了书本上的矩阵相乘的例子,这样开始写代码,然后很自然地出现了各种问题。

以下的内容供大家学习参考,有问题可以留言与我反馈。

开始学着使用 CUDA,实现一个矩阵乘法运算。

首先我们要定义一个矩阵的结构体,话说CUDA是否支持结构体作为设备端的函数的参数呢?

不妨都一股脑试验一下。

1、安装CUDA 5.0

在NVIDIA CUDA官方网站下载对应自己操作系统的最新版 CUDA,

为什么要最新版本呢,一方面技术一直更新,自己要保持潮流嘛~(程序员的苦逼之处也在此)

更新的版本一般实现功能上更加完善,BUG更少。

废话不多说,去这里下载:

最新版的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;
}

posted on 2015-07-12 12:08  袁军峰  阅读(447)  评论(0编辑  收藏  举报

导航