CUDA编程记之二 详解GPU中各种内存属性及实战演练

系列文章目录

CUDA编程记之一基本使用及线程、同步、存储器

---

---

前言

像之前的文章，也只能说讲了一下简单应用，其实离实际应用还有很大距离，这篇再细讲讲存储器和应用示例

一、存储器和内存

图中所示为GPU中的存储结构，L1 、 L2为缓存区域，全局内存访问很慢，但所有位置都可以访问。共享内存访问比全局内存快100倍左右，但是是块内可访问，不同块之间的共享内存不同。
本地内存访问也很快，但是只有64kb左右。溢出时候会占用寄存器内存，也是很慢的。本地内存一般是指函数内声明的局部变量。

所以理论上来说，对全局内存的操作越少越好，尽量都在共享内存中完成。具体在后面示例中就可以学习了

二、矩阵点积

#include "stdio.h"
#include<iostream>
#include <cuda.h>
#include <cuda_runtime.h>
#define N 1024
#define threadsPerBlock 512

/*
演示了怎么使用共享内存，以及利用索引累加，总而言之能并行就并行
*/

__global__ void gpu_dot(float *d_a, float *d_b, float *d_c) {
	// 定义与块内线程数相同的共享内存数组
	__shared__ float partial_sum[threadsPerBlock];
	// 定义一个总的索引
	int tid = threadIdx.x + blockIdx.x * blockDim.x;
	// 计算块内索引
	int index = threadIdx.x;
	//Calculate Partial Sum
	float sum = 0;
	// 通过这种方法，当线程数x块数小于要运算的总次数时候，也可以计算
	while (tid < N)
	{
		// 因为是点积，所以是对应元素相乘
		sum += d_a[tid] * d_b[tid];
		// 每次累加总线程数
		tid += blockDim.x * gridDim.x;
	}

	// 把每个线程的结果赋值到共享内存
	partial_sum[index] = sum;

	/*
	相当于每个矩阵有n个元素，你能开启的就是blockDim.x * gridDim.x 个线程，所以根据总线程数取模，部分线程会运算多几次，但是肯定不会有元素重复运算，接着每个线程全都把结果赋予到共享内存，下一步按道理是求总结果
	那是不是还需要共享内存数组，就一个共享内存总数不行吗？
	*/


	// 同步一次，保证到此为止，所以线程都计算并将结果存到共享内存
	__syncthreads();

	// 还是不确定和原子操作哪个更快，现在看起来是这种方式累加比较快

	// 累加到只剩一个
	int i = blockDim.x / 2;
	while (i != 0) {
		if (index < i)
			partial_sum[index] += partial_sum[index + i];
		__syncthreads();
		i /= 2;
	}
	// 把每个块内的结果赋予到d_c
	// 原书的意思是这里可以使用原子操作直接累加结束，然后就传回一个数
	if (index == 0)
		d_c[blockIdx.x] = partial_sum[0];
}


int main(void) {
	// 定义主机指针，应该是partial_sum用于传出，h_c表示结果
	float *h_a, *h_b, h_c, *partial_sum;
	// 定义设备指针
	float *d_a, *d_b, *d_partial_sum;
	//取整，此时block_calc * threadsPerBlock 总数大于N
	int block_calc = (N + threadsPerBlock - 1) / threadsPerBlock;
	// 最大取到32个块
	int blocksPerGrid = (32 < block_calc ? 32 : block_calc);
	// 主机指针也申请内存
	h_a = (float*)malloc(N * sizeof(float));
	h_b = (float*)malloc(N * sizeof(float));
	partial_sum = (float*)malloc(blocksPerGrid * sizeof(float));

	// 设备指针申请内存
	cudaMalloc((void**)&d_a, N * sizeof(float));
	cudaMalloc((void**)&d_b, N * sizeof(float));
	cudaMalloc((void**)&d_partial_sum, blocksPerGrid * sizeof(float));

	// A矩阵按0到N-1递增，B矩阵全都是2
	// 那预估结果就是(0+N-1)/2*N*2
	for (int i = 0; i<N; i++) {
		h_a[i] = i;
		h_b[i] = 2;
	}

	//  把主机的值拷贝到设备
	cudaMemcpy(d_a, h_a, N * sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(d_b, h_b, N * sizeof(float), cudaMemcpyHostToDevice);
	// 运算过程
	gpu_dot << <blocksPerGrid, threadsPerBlock >> >(d_a, d_b, d_partial_sum);

	// 把值拷贝会主机
	cudaMemcpy(partial_sum, d_partial_sum, blocksPerGrid * sizeof(float), cudaMemcpyDeviceToHost);

	// 总结果置0
	h_c = 0;
	// 累加
	for (int i = 0; i<blocksPerGrid; i++) {
		h_c += partial_sum[i];
	}
	printf("The computed dot product is: %f\n", h_c);
	// 这边是比对下结果
#define cpu_sum(x) (x*(x+1))

	// 浮点数不能用==，还是检测绝对值之差是否合格比较ok
	if (h_c == cpu_sum((float)(N - 1)))
	{
		printf("The dot product computed by GPU is correct\n");
	}
	else
	{
		printf("Error in dot product computation");
	}


	// 该释放还是释放
	cudaFree(d_a);
	cudaFree(d_b);
	cudaFree(d_partial_sum);
	free(h_a);
	free(h_b);
	free(partial_sum);

	system("pause");
}

其实比较好理解，点积就是两个矩阵对应元素之积，那么每个线程就是相乘即可，如果能同时开启的线程不够，那么可以通过while循环来实现倍率的运算。
原书的意思应该是运算完，结果怎么处理，是传回一个数组，由cpu累加还是直接在GPU全部完成。
ps：块内的结果累加使用代码中的方法应该是由于原子操作直接累加的，要后续才能验证

矩阵乘法

// 矩阵乘法
#include "stdio.h"
#include<iostream>
#include <cuda.h>
#include <cuda_runtime.h>
#include <math.h>
#define TILE_SIZE 2


// 不使用共享内存的菜鸡核函数
__global__ void gpu_Matrix_Mul_nonshared(float *d_a, float *d_b, float *d_c, const int size)
{
	int row, col;
	col = TILE_SIZE * blockIdx.x + threadIdx.x;
	row = TILE_SIZE * blockIdx.y + threadIdx.y;

	for (int k = 0; k< size; k++)
	{
		// 相比下面的内容，这里就实在是太简单了
		d_c[row*size + col] += d_a[row * size + k] * d_b[k * size + col];
	}
}

// 使用共享内存的核函数组长（不同于挥舞皮鞭的督工，是确实有更厉害）
__global__ void gpu_Matrix_Mul_shared(float *d_a, float *d_b, float *d_c, const int size)
{
	int row, col;
	// 定义共享二维数组
	__shared__ float shared_a[TILE_SIZE][TILE_SIZE];
	__shared__ float shared_b[TILE_SIZE][TILE_SIZE];
	// 第一遍看都没看懂，呆了十来分钟
	// 因为是二维的块和线程，使用可以理解为把块和线程组合好以后能画出一个和矩阵同形状的网格
	// 所以确实是col和row，只不过第一次确实要想一下，我也不好说是我蠢
	col = TILE_SIZE * blockIdx.x + threadIdx.x;
	row = TILE_SIZE * blockIdx.y + threadIdx.y;

	// 这里又卡了我一个小时，理解为什么要这么操作，感觉相对来说自己的思想不够高维
	// 第一层for循环遍历所有block，最好先理解一个block内发生了什么
	// 最后我通过直接写出几个格子的内容，确认是正确的，但还是奇怪他们是怎么想到这么写的
	// 如果是我，我会怎么写
	for (int i = 0; i< size / TILE_SIZE; i++)
	{
		shared_a[threadIdx.y][threadIdx.x] = d_a[row* size + (i*TILE_SIZE + threadIdx.x)];
		shared_b[threadIdx.y][threadIdx.x] = d_b[(i*TILE_SIZE + threadIdx.y) * size + col];
		__syncthreads();
		for (int j = 0; j<TILE_SIZE; j++)
			d_c[row*size + col] += shared_a[threadIdx.y][j] * shared_b[j][threadIdx.x];
		__syncthreads();

	}

	/*
	关于自己会怎么写，我觉得首先最简单的就是直接定义三个共享内存的二维数组，按照之前那个非共享的方式进行运算
	即只有一个block，内部有个二维的线程组，然后运算全都在共享内存上进行，最后再赋值

	下面是改进计划：
	1.不能只有一个block
	2.如果是奇数的边长怎么办

	他的方式看起来很简洁，但感觉也不好说是最好

	如果还是按他的这么多个BLOCK,这样的THREAD结构，那我怎么样让共享内存生效

	叮叮叮，顿悟了

	可以这么理解，把矩阵分解成n个block，为了便于理解，还是使用二维block，然后二维thread
	这样组合起来还是和矩阵形式相同的网格，非常好理解

	这时你要通过共享内存来加速，还是按这里的2x2小矩阵看，求某个点的结果值，是A矩阵的整行和B矩阵的整列的点积

	说明2x2小矩阵中的同行元素共用a的同行，同列元素共用b的同列。
	可以理解为在原来的AB矩阵上套了这个2x2的范围，然后为了运算完所有需要的值，窗口在A上横向滑动，在B上纵向滑动，
	内部每次计算累加
	总之就是还比较精妙，不得不感慨，我的数学功底比牛顿 爱因斯坦应该还是要差一点
	感觉要是要写一个加速的不定尺寸的矩阵乘法已经很复杂了。

	*/
}

// main routine
int main()
{
	// 矩阵的大小
	const int size = 4;
	// 定义主机数组，4x4方阵
	float h_a[size][size], h_b[size][size], h_result[size][size];
	// 定义设备指针
	float *d_a, *d_b, *d_result;
	// A每个元素等于自己的行数,B每个元素都等于自己的列
	for (int i = 0; i<size; i++)
	{
		for (int j = 0; j<size; j++)
		{
			h_a[i][j] = i;
			h_b[i][j] = j;
		}
	}

	// 为设备指针申请内存
	cudaMalloc((void **)&d_a, size*size * sizeof(int));
	cudaMalloc((void **)&d_b, size*size * sizeof(int));
	cudaMalloc((void **)&d_result, size*size * sizeof(int));


	// 把主机上的值拷贝到设备

	cudaMemcpy(d_a, h_a, size*size * sizeof(int), cudaMemcpyHostToDevice);
	cudaMemcpy(d_b, h_b, size*size * sizeof(int), cudaMemcpyHostToDevice);

	// 怎么看起来要使用三维的块，总共4个块
	dim3 dimGrid(size / TILE_SIZE, size / TILE_SIZE, 1);
	// 每块有4个线程
	dim3 dimBlock(TILE_SIZE, TILE_SIZE, 1);
	//gpu_Matrix_Mul_nonshared << <dimGrid, dimBlock >> > (d_a, d_b, d_result, size);

	gpu_Matrix_Mul_shared << <dimGrid, dimBlock >> > (d_a, d_b, d_result, size);

	cudaMemcpy(h_result, d_result, size*size * sizeof(int), cudaMemcpyDeviceToHost);
	printf("The result of Matrix multiplication is: \n");

	for (int i = 0; i< size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			printf("%f   ", h_result[i][j]);
		}
		printf("\n");
	}
	cudaFree(d_a);
	cudaFree(d_b);
	cudaFree(d_result);
	return 0;
}

---

总结

提示：这里对文章进行总结：
总之有点难，我今天要消化一下