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

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

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

开始学着使用 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语言基础的人能够轻松看懂这个函数而不需要我的注释。

// 矩阵乘法，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的定义如下：

// 矩阵结构体
typedef struct MartixTag
{
	float *ele;	// 一维数组，保存矩阵的元素
	size_t col;	// 矩阵的列数
	size_t row;	// 矩阵的行数
	size_t stride;	// 矩阵一行数据的数量，为了方便内存对齐，2^n
}Matrix;

实际测试表明，上述代码中，采用指针传递结构体参数时，运算结果会不正确。

这点在我们以后写CUDA函数的时候应该注意，我的理解是这样子，不知道大家怎么看：

// （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、附录代码

最后附上我的一个成果，按照书上的分块计算矩阵乘法代码，经本人测试无误：

__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;
}

调用方式如下（对相关参数初始化完毕后执行）： 

// 启用 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源代码计算矩阵乘积，也分享给大家：

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;
}