CUDA编程的理论部分可以参考模型部署篇 中的GPU 的 CUDA 编程方法。
虽然CUDA有很多的C代码，这里我们主要以C++为主。一个完整的CUDA程序，需要经历7个步骤
- 设置显卡设备
- 分配显存空间
- 从内存到显存拷贝数据
- 执行CUDA并行函数
- CUDA函数结束后，将结果从显存拷贝回内存
- 释放显存空间
- 设备重置
如果是单GPU的话可以省略1跟7两个步骤。
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime_api.h>
#include <iostream>
using namespace std;
/* 核函数 */
__global__ void kernelFunc(float *a) {
a[threadIdx.x] = 1;
}
int main(int argc, char **argv) {
//获取GPU的数量
int gpuCount = -1;
cudaGetDeviceCount(&gpuCount);
cout << gpuCount << endl;
if (gpuCount < 0) {
cout << "no device" << endl;
}
//设置显卡设备
cudaSetDevice(0);
//分配显存空间
float *aGpu;
cudaMalloc((void**)&aGpu, 16 * sizeof(float));
//从内存到显存拷贝数据
float a[16] = {0};
cudaMemcpy(aGpu, a, 16 * sizeof(float), cudaMemcpyHostToDevice);
//执行CUDA并行函数
kernelFunc <<<1, 16>> >(aGpu);
//CUDA函数结束后，将结果从显存拷贝回内存
cudaMemcpy(a, aGpu, 16 * sizeof(float), cudaMemcpyDeviceToHost);
for (int i = 0; i < 16; i++) {
printf("%f", a[i]);
}
printf("\n");
//释放显存空间
cudaFree(aGpu);
//设备重置
cudaDeviceReset();
//获取设备的属性
cudaDeviceProp prop;
cudaGetDeviceProperties(&prop, 0);
//块最大线程数
printf("maxThreadsPerBlock: %d\n", prop.maxThreadsPerBlock);
//块维度最大值
printf("maxThreadsDim: %d\n", prop.maxThreadsDim[0]);
//Grid各维度最大值
printf("maxGridSize: %d\n", prop.maxGridSize[0]);
//常量内存的大小
printf("totalConstMem: %d\n", prop.totalConstMem);
//时钟频率
printf("clockRate: %d\n", prop.clockRate);
//GPU是否集成
printf("initegrated: %d\n", prop.integrated);
return 0;
}
CUDA的源码文件以.cu为后缀，编译命令如下(需要先安装CUDA，安装方式可以参考乌班图安装Pytorch、Tensorflow Cuda环境 )
nvcc main.cu -o main
运行结果
1
1.0000001.0000001.0000001.0000001.0000001.0000001.0000001.0000001.0000001.0000001.0000001.0000001.0000001.0000001.0000001.000000
maxThreadsPerBlock: 1024
maxThreadsDim: 1024
maxGridSize: 2147483647
totalConstMem: 65536
clockRate: 921600
initegrated: 1
基本运算
加法运算
#include <iostream>
using namespace std;
//核函数
__global__ void add(int* a, int* b, int* c, int num) {
//获取GPU线程的Id
int i = threadIdx.x;
if (i < num) {
c[i] = a[i] + b[i];
}
}
int main() {
int num = 10;
int a[num], b[num], c[num];
int *a_g, *b_g, *c_g;
for (int i = 0; i < num; i++) {
a[i] = i;
b[i] = i * i;
}
//在GPU上分配空间
cudaMalloc((void**)&a_g, num * sizeof(int));
cudaMalloc((void**)&b_g, num * sizeof(int));
cudaMalloc((void**)&c_g, num * sizeof(int));
//从CPU上拷贝数据到GPU上
cudaMemcpy(a_g, a, num * sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(b_g, b, num * sizeof(int), cudaMemcpyHostToDevice);
add<<<1, num>>>(a_g, b_g, c_g, num);
//从GPU上拷贝到CPU上
cudaMemcpy(c, c_g, num * sizeof(int), cudaMemcpyDeviceToHost);
for (int i = 0; i < num; i++) {
cout << a[i] << "+" << b[i] << "=" << c[i] << endl;
}
return 0;
}
编译
nvcc main.cu -o main
运行结果
0+0=0
1+1=2
2+4=6
3+9=12
4+16=20
5+25=30
6+36=42
7+49=56
8+64=72
9+81=90
用GPU写卷积核
#include <iostream>
using namespace std;
//获取GPU的线程数
int getThreadNum() {
cudaDeviceProp prop;
//GPU数
int count;
cudaGetDeviceCount(&count);
cout << "gpu num: " << count << endl;
//获取GPU的参数
cudaGetDeviceProperties(&prop, 0);
//每个线程块的最大线程数
cout << "max thread num: " << prop.maxThreadsPerBlock << endl;
//查看三个网格各有多少个线程块
cout << "grid dimensions:" << prop.maxGridSize[0] << "," << prop.maxGridSize[1] << "," << prop.maxGridSize[2] << endl;
return prop.maxThreadsPerBlock;
}
//卷积运算
__global__ void conv(float *img, float *kernel, float *result, int width, int height, int kernelSize) {
//当前线程值
int ti = threadIdx.x;
//当前块值
int bi = blockIdx.x;
int id = (bi * blockDim.x + ti);
if (id >= width * height) {
return;
}
//中心点
int row = id / width;
int col = id % width;
for(int i = 0; i < kernelSize; i++) {
for(int j = 0; j < kernelSize; j++) {
float imgValue = 0;
int curRow = row - kernelSize / 2 + i;
int curCol = col - kernelSize / 2 + j;
if (curRow < 0 || curCol < 0 || curRow >= height || curCol >= width) {
} else {
imgValue = img[curRow * width + curCol];
}
result[id] += kernel[i * kernelSize + j] * imgValue;
}
}
}
int main() {
int width = 1920;
int height = 1080;
//创建一张图片
float *img = new float[width * height];
float *result = new float[width * height];
for(int row = 0; row < height; row++) {
for(int col = 0; col < width; col++) {
img[col + row * width] = (col + row) % 256;
}
}
//卷积核的大小
int kernelSize = 3;
//卷积核
float *kernel = new float[kernelSize * kernelSize];
for(int i = 0; i < kernelSize * kernelSize; i++) {
kernel[i] = i % kernelSize - 1;
}
//创建GPU上的对象
float *img_g;
float *kernel_g;
float *result_g;
//在GPU上分配空间
cudaMalloc((void**)&img_g, width * height * sizeof(float));
cudaMalloc((void**)&kernel_g, kernelSize * kernelSize * sizeof(float));
cudaMalloc((void**)&result_g, width * height * sizeof(float));
//从CPU上将数据拷贝到GPU上
cudaMemcpy(img_g, img, width * height * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(kernel_g, kernel, kernelSize * kernelSize * sizeof(float), cudaMemcpyHostToDevice);
int threadNum = getThreadNum();
int blockNum = (width * height - 0.5) / threadNum + 1;
conv<<<blockNum, threadNum>>>(img_g, kernel_g, result_g, width, height, kernelSize);
cudaMemcpy(result, result_g, width * height * sizeof(float), cudaMemcpyDeviceToHost);
cout << "img" << endl;
//打印图图片的10*10的值
for(int row = 0; row < 10; row++) {
for(int col = 0; col < 10; col++) {
cout << img[col + row * width] << " ";
}
cout << endl;
}
cout << "kernel" << endl;
//打印卷积核
for(int row = 0; row < kernelSize; row++) {
for(int col = 0; col < kernelSize; col++) {
cout << kernel[col + row * kernelSize] << " ";
}
cout << endl;
}
cout << "result" << endl;
for(int row = 0; row < 10; row++) {
for(int col = 0; col < 10; col++) {
cout << result[col + row * width] << " ";
}
cout << endl;
}
return 0;
}
运行结果
gpu num: 1
max thread num: 1024
grid dimensions:2147483647,65535,65535
img
0 1 2 3 4 5 6 7 8 9
1 2 3 4 5 6 7 8 9 10
2 3 4 5 6 7 8 9 10 11
3 4 5 6 7 8 9 10 11 12
4 5 6 7 8 9 10 11 12 13
5 6 7 8 9 10 11 12 13 14
6 7 8 9 10 11 12 13 14 15
7 8 9 10 11 12 13 14 15 16
8 9 10 11 12 13 14 15 16 17
9 10 11 12 13 14 15 16 17 18
kernel
-1 0 1
-1 0 1
-1 0 1
result
3 4 4 4 4 4 4 4 4 4
6 6 6 6 6 6 6 6 6 6
9 6 6 6 6 6 6 6 6 6
12 6 6 6 6 6 6 6 6 6
15 6 6 6 6 6 6 6 6 6
18 6 6 6 6 6 6 6 6 6
21 6 6 6 6 6 6 6 6 6
24 6 6 6 6 6 6 6 6 6
27 6 6 6 6 6 6 6 6 6
30 6 6 6 6 6 6 6 6 6
- 网格(Grid)、线程块(Block)、线程(Thread)
CUDA的软件架构由网格（Grid）、线程块（Block）和线程（Thread）组成，相当于把GPU上的计算单元分为若干（2~3）个网格，每个网格内包含若干（65535）个线程块，每个线程块包含若干（1024）个线程。
- 线程索引的计算公式
我们写核函数的时候第一步是一定要得到线程索引的，但是线程索引会根据情况不同而不同，并不一定就是threadIdx.x
一个Grid可以包含多个Blocks，Blocks的组织方式可以是一维的，二维或者三维的。block包含多个Threads，这些Threads的组织方式也可以是一维，二维或者三维的。
CUDA中每一个线程都有一个唯一的标识ID—ThreadIdx，这个ID随着Grid和Block的划分方式的不同而变化，这里给出Grid和Block不同划分方式下线程索引ID的计算公式。
1、grid划分成1维，block划分成1维
int threadId = blockIdx.x * blockDim.x + threadIdx.x;
在我们卷积核的案例中就是使用的该方法来获取的线程索引。但是在加法运算中则直接使用的是threadIdx.x。这是因为在加法运算中我们的数据量很小，算的是两个长度为10的数组的加法。我们知道一个线程块有1024个线程，完全够用来计算20个数的相加，所以这里只需要使用一个线程块就够了。而在卷积核的计算中，我们要计算的是一个1920*1080=2073600的图像中所有像素的值，2073600已经远远超过了1024，一个线程块肯定是不够用的，故而需要使用多个(2073600/1024=2025)线程块来进行处理。而一个网格内有65535个线程块，2025<65535，故而在一个网格内就可以计算。故上式我们就可以理解成第几个线程块的第几个线程合起来组成了我们需要的线程索引。
在上面的两个例程中会有这样两段代码
add<<<1, num>>>(a_g, b_g, c_g, num);
conv<<<blockNum, threadNum>>>(img_g, kernel_g, result_g, width, height, kernelSize);
这里<<<a,b>>>中间的两个数a,b就代表我们运行该核函数需要开多少个线程块和线程数。
2、grid划分成1维，block划分成2维
int threadId = blockIdx.x * blockDim.x * blockDim.y+ threadIdx.y * blockDim.x + threadIdx.x;
3、grid划分成1维，block划分成3维
int threadId = blockIdx.x * blockDim.x * blockDim.y * blockDim.z
+ threadIdx.z * blockDim.y * blockDim.x
+ threadIdx.y * blockDim.x + threadIdx.x;
4、grid划分成2维，block划分成1维
int blockId = blockIdx.y * gridDim.x + blockIdx.x;
int threadId = blockId * blockDim.x + threadIdx.x;
5、grid划分成2维，block划分成2维
int blockId = blockIdx.x + blockIdx.y * gridDim.x;
int threadId = blockId * (blockDim.x * blockDim.y)
+ (threadIdx.y * blockDim.x) + threadIdx.x;
6、grid划分成2维，block划分成3维
int blockId = blockIdx.x + blockIdx.y * gridDim.x;
int threadId = blockId * (blockDim.x * blockDim.y * blockDim.z)
+ (threadIdx.z * (blockDim.x * blockDim.y))
+ (threadIdx.y * blockDim.x) + threadIdx.x;
7、grid划分成3维，block划分成1维
int blockId = blockIdx.x + blockIdx.y * gridDim.x
+ gridDim.x * gridDim.y * blockIdx.z;
int threadId = blockId * blockDim.x + threadIdx.x;
8、grid划分成3维，block划分成2维
int blockId = blockIdx.x + blockIdx.y * gridDim.x
+ gridDim.x * gridDim.y * blockIdx.z;
int threadId = blockId * (blockDim.x * blockDim.y)
+ (threadIdx.y * blockDim.x) + threadIdx.x;
9、grid划分成3维，block划分成3维
int blockId = blockIdx.x + blockIdx.y * gridDim.x
+ gridDim.x * gridDim.y * blockIdx.z;
int threadId = blockId * (blockDim.x * blockDim.y * blockDim.z)
+ (threadIdx.z * (blockDim.x * blockDim.y))
+ (threadIdx.y * blockDim.x) + threadIdx.x;
数组求和
#include <iostream>
using namespace std;
__global__ void sum(float *a, float *result) {
int tid = threadIdx.x;
//定义一块同一个线程块的共享内存
__shared__ float sData[16];
//将每一个线程的值都放入到共享内存中
sData[tid] = a[tid];
//同步所有的线程全部写入共享内存
__syncthreads();
for(int i = 8; i > 0; i /= 2) {
if (tid < i) {
sData[tid] += sData[tid + i];
}
__syncthreads();
}
if (tid == 0) {
result[0] = sData[0];
}
}
int main() {
float a[16];
float result[1];
for(int i = 0; i < 16; i++) {
a[i] = i * (i + 1);
}
float *a_g;
cudaMalloc((void**)&a_g, 16 * sizeof(float));
cudaMemcpy(a_g, a, 16 * sizeof(float), cudaMemcpyHostToDevice);
float *result_g;
cudaMalloc((void**)&result_g, 1 * sizeof(float));
sum<<<1,16>>>(a_g, result_g);
cudaMemcpy(result, result_g, 1 * sizeof(float), cudaMemcpyDeviceToHost);
cout << result[0] << endl;
return 0;
}
运行结果
1360
- 共享内存
GPU 包括板载和片上两种内存：全局内存是较大的板载内存，具有相对较高的延迟；共享内存是较小的片上内存，具有相对较低的延迟，可提供比全局内存更高的带宽，可当作一个可编程管理的缓存。
共享内存的应用场景：
- 块内线程通信的通道；
- 用于全局内存数据的可编程管理的缓存；
- 高速暂存存储器，用于转换数据以优化全局内存访问模式
共享内存（shared memory， SMEM） 是GPU的一个关键部件。物理上，每个SM都有 一个小的低延迟内存池， 这个内存池被当前正在该SM上执行的线程块中的所有线程所共享。SMEM使同一个线程块中的线程能够互相协作，便于片上数据重用，降低核函数所需的全局内存带宽。
全局内存的所有加载和存储请求都要经过二级缓存， 这是SM单元之间数据统一的基本点。 相较于二级缓存和全局内存，共享内存和一级缓存在物理上更接近SM。
当每个线程块开始执行时，会分配给它一定数量的共享内存。 这个共享内存的地址空间被线程块中所有的线程共享。它和创建时所在的线程块具有相同生命周期。线程束访问共享内存：最好，访问请求在1个事务中完成；最坏，在32个不同的事务中顺序执行。
共享内存被SM中的所有常驻线程块划分，是设备并行性的关键资源。一个核函数使用的共享内存越多，处于并发活跃状态的线程块就越少，因为虽然一个SM可以有多个活跃线程块，如果核函数占用的共享内存过多，就挤占了其线程块的资源，没有分配到共享内存资源的线程块就只有等待，从而减少了活跃线程块的数量。
CUDA中允许手动管理共享内存，在数据布局上提供更多的细粒度控制，简化了代码优化。
共享内存变量用__shared__修饰符声明，可以静态或动态地分配内存；可以本地或全局声明的CUDA核函数。CUDA支持1维、2维和3维共享内存数组的声明。
GPU与CPU的速度比较
- IO频繁简单操作
#include <iostream>
#include <time.h>
#include <sys/time.h>
using namespace std;
__global__ void sum(float *a, float *result) {
int tid = threadIdx.x;
//定义一块同一个线程块的共享内存
__shared__ float sData[16];
//将每一个线程的值都放入到共享内存中
sData[tid] = a[tid];
//同步所有的线程全部写入共享内存
__syncthreads();
for(int i = 8; i > 0; i /= 2) {
if (tid < i) {
sData[tid] += sData[tid + i];
}
__syncthreads();
}
if (tid == 0) {
result[0] = sData[0];
}
}
void sumCpu(float *a, float *result) {
result[0] = 0;
for(int i = 0; i < 16; i++) {
result[0] += a[i];
}
}
int main() {
float a[16];
float result[1];
for(int i = 0; i < 16; i++) {
a[i] = i * (i + 1);
}
float *a_g;
cudaMalloc((void**)&a_g, 16 * sizeof(float));
cudaMemcpy(a_g, a, 16 * sizeof(float), cudaMemcpyHostToDevice);
float *result_g;
cudaMalloc((void**)&result_g, 1 * sizeof(float));
struct timeval startTime, endTime;
int sunM = 1000000;
gettimeofday(&startTime, NULL);
for(int i = 0; i < sunM; i++) {
sum<<<1,16>>>(a_g, result_g);
}
gettimeofday(&endTime, NULL);
cout << "cuda use time: " << (endTime.tv_sec - startTime.tv_sec) * 1000000 + (endTime.tv_sec - startTime.tv_sec) << endl;
cudaMemcpy(result, result_g, 1 * sizeof(float), cudaMemcpyDeviceToHost);
cout << result[0] << endl;
gettimeofday(&startTime, NULL);
for(int i = 0; i< sunM; i++) {
sumCpu(a, result);
}
gettimeofday(&endTime, NULL);
cout << "cpu use time: " << (endTime.tv_sec - startTime.tv_sec) * 1000000 + (endTime.tv_sec - startTime.tv_sec) << endl;
cout << result[0] << endl;
return 0;
}
运行结果
cuda use time: 2000002
1360
cpu use time: 0
1360
在这种情况下，我们可以看出GPU的运算速度要远远慢于CPU的，因为这是一个16个数相加的计算，运行1000000次，对于GPU来说有着大量的从CPU拷贝到GPU的数据操作以及从GPU拷贝回CPU的数据操作，这个过程是比较耗时的，也就是时间都浪费在了IO操作上。而对于CPU运算则没有这种大量的IO操作，所以CPU的运算速度要快的多。故而GPU并不适合于大量频繁的IO操作的简单计算。
- 低频大数据量计算
#include <iostream>
#include <time.h>
#include <sys/time.h>
using namespace std;
__global__ void add(int *a, int *b, int *result, int num) {
int i;
int tid = threadIdx.x;
int bid = blockIdx.x;
int threadNum = blockDim.x;
i = bid * threadNum + tid;
if (i < num) {
result[i] = a[i] + b[i];
}
}
void addCpu(int *a, int *b, int *result, int num) {
for(int i = 0; i < num; i++) {
result[i] = a[i] + b[i];
}
}
int main() {
int num = 512000;
int threadNum = 128;
int blockNum = 40;
int a[num], b[num], result[num];
for(int i = 0; i < num; i++) {
a[i] = i;
b[i] = i * i;
}
int *a_g, *b_g, *result_g;
cudaMalloc((void**)&a_g, num * sizeof(int));
cudaMalloc((void**)&b_g, num * sizeof(int));
cudaMalloc((void**)&result_g, num * sizeof(int));
cudaMemcpy(a_g, a, num * sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(b_g, b, num * sizeof(int), cudaMemcpyHostToDevice);
struct timeval startTime, endTime;
int loopNum = 10000;
gettimeofday(&startTime, NULL);
for(int i = 0; i < loopNum; i++) {
add<<<blockNum,threadNum>>>(a_g, b_g, result_g, num);
}
gettimeofday(&endTime, NULL);
cout << "cuda use time: " << (endTime.tv_sec - startTime.tv_sec) * 1000000 + (endTime.tv_sec - startTime.tv_sec) << endl;
cudaMemcpy(result, result_g, num * sizeof(int), cudaMemcpyDeviceToHost);
for(int i = 0; i < 50; i++) {
cout << result[i] << " ";
}
cout << endl;
gettimeofday(&startTime, NULL);
for(int i = 0; i< loopNum; i++) {
addCpu(a, b, result, num);
}
gettimeofday(&endTime, NULL);
cout << "cpu use time: " << (endTime.tv_sec - startTime.tv_sec) * 1000000 + (endTime.tv_sec - startTime.tv_sec) << endl;
for(int i = 0; i < 50; i++) {
cout << result[i] << " ";
}
cout << endl;
return 0;
}
运行结果
cuda use time: 0
0 2 6 12 20 30 42 56 72 90 110 132 156 182 210 240 272 306 342 380 420 462 506 552 600 650 702 756 812 870 930 992 1056 1122 1190 1260 1332 1406 1482 1560 1640 1722 1806 1892 1980 2070 2162 2256 2352 2450
cpu use time: 4000004
0 2 6 12 20 30 42 56 72 90 110 132 156 182 210 240 272 306 342 380 420 462 506 552 600 650 702 756 812 870 930 992 1056 1122 1190 1260 1332 1406 1482 1560 1640 1722 1806 1892 1980 2070 2162 2256 2352 2450
这里是每次计算2个512000长度的数组相加，执行10000次，10000次相对之前的1000000次要小了很多，但是单次的计算量要大了很多。这里我们可以看到GPU的运行时间要远远小于CPU。
原子操作
我们把之前的数组求和修改一下
#include <iostream>
using namespace std;
__global__ void sum(float *a, float *result) {
int tid = threadIdx.x;
result[0] = 0;
__syncthreads();
result[0] += a[tid];
}
int main() {
float a[16];
float result[1];
for(int i = 0; i < 16; i++) {
a[i] = i * (i + 1);
}
float *a_g;
cudaMalloc((void**)&a_g, 16 * sizeof(float));
cudaMemcpy(a_g, a, 16 * sizeof(float), cudaMemcpyHostToDevice);
float *result_g;
cudaMalloc((void**)&result_g, 1 * sizeof(float));
sum<<<1,16>>>(a_g, result_g);
cudaMemcpy(result, result_g, 1 * sizeof(float), cudaMemcpyDeviceToHost);
cout << result[0] << endl;
return 0;
}
在这里我们是把每一个线程索引的a值都加一遍，但是执行结果却是
0
这明显跟我们之前算得的1360不符。究其原因就是它这里是线程不安全的计算，对于多个线程同时操作一块内存空间的时候，线程安全问题是由其需要注意的。
- 计算直方图
#include <iostream>
using namespace std;
__global__ void get_hist(float *a, int *hist) {
int tid = threadIdx.x;
int bid = blockIdx.x;
int idx = tid + bid * blockDim.x;
//原子+1操作
atomicAdd(&hist[(int)a[idx]], 1);
}
int main() {
int size = 320000;
float *a = new float[size];
int length = 255;
for(int i = 0; i < size; i++) {
a[i] = i * (i + 1) % length;
}
int hist[length] = {0};
float *a_g;
int *hist_g;
cudaMalloc((void**)&a_g, size * sizeof(float));
cudaMalloc((void**)&hist_g, length * sizeof(int));
cudaMemcpy(a_g, a, size * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(hist_g, hist, length * sizeof(int), cudaMemcpyHostToDevice);
get_hist<<<size / 512, 512>>>(a_g, hist_g);
cudaMemcpy(hist, hist_g, length * sizeof(int), cudaMemcpyDeviceToHost);
cout << "gpu" << endl;
for(int i = 0; i < length; i++) {
cout << hist[i] << " ";
}
cout << endl;
int hist1[length] = {0};
for(int i = 0; i < size; i++) {
hist1[(int)a[i]] += 1;
}
cout << "cpu" << endl;
for(int i = 0; i < length; i++) {
cout << hist1[i] << " ";
}
cout << endl;
return 0;
}
运行结果
gpu
2855 769 1312 642 747 1322 1110 519 551 600 669 807 1820 443 476 522 606 1134 449 635 1174 709 362 377 454 481 513 556 681 812 1822 387 408 467 620 696 1215 515 552 600 669 805 1818 387 415 459 536 1055 383 478 557 1076 375 408 462 552 921 1605 211 218 226 233 241 250 370 386 405 530 568 609 682 1003 1151 447 467 490 524 572 629 842 1419 1213 625 665 778 860 1016 2099 716 848 1849 410 426 448 467 487 519 554 596 711 798 956 2053 676 825 1856 478 994 330 498 1002 285 297 303 317 332 347 396 474 504 553 634 1162 478 550 1066 362 392 425 485 573 948 1637 242 259 271 291 313 378 440 754 891 274 284 291 302 314 383 400 420 441 472 513 592 695 1263 1049 343 365 452 507 788 463 531 661 1665 229 249 297 470 974 259 380 399 422 466 731 434 455 484 531 603 1108 428 553 579 609 751 808 927 1032 1312 2173 846 985 2002 576 612 690 772 1296 601 698 751 878 1002 1404 2087 711 863 1895 517 1030 330 367 437 945 337 363 416 464 733 404 503 585 1117 455 960 266 279 290 301 422 444 518 550 587 692 821 963 1992 582 857 522 677 754 1271 572 611 705 771 1060 813 930 1085
cpu
2855 769 1312 642 747 1322 1110 519 551 600 669 807 1820 443 476 522 606 1134 449 635 1174 709 362 377 454 481 513 556 681 812 1822 387 408 467 620 696 1215 515 552 600 669 805 1818 387 415 459 536 1055 383 478 557 1076 375 408 462 552 921 1605 211 218 226 233 241 250 370 386 405 530 568 609 682 1003 1151 447 467 490 524 572 629 842 1419 1213 625 665 778 860 1016 2099 716 848 1849 410 426 448 467 487 519 554 596 711 798 956 2053 676 825 1856 478 994 330 498 1002 285 297 303 317 332 347 396 474 504 553 634 1162 478 550 1066 362 392 425 485 573 948 1637 242 259 271 291 313 378 440 754 891 274 284 291 302 314 383 400 420 441 472 513 592 695 1263 1049 343 365 452 507 788 463 531 661 1665 229 249 297 470 974 259 380 399 422 466 731 434 455 484 531 603 1108 428 553 579 609 751 808 927 1032 1312 2173 846 985 2002 576 612 690 772 1296 601 698 751 878 1002 1404 2087 711 863 1895 517 1030 330 367 437 945 337 363 416 464 733 404 503 585 1117 455 960 266 279 290 301 422 444 518 550 587 692 821 963 1992 582 857 522 677 754 1271 572 611 705 771 1060 813 930 1085
这里的atomicAdd(&hist[(int)a[idx]], 1);保证了多个线程无法对hist同时进行加1操作导致结果错误，最终GPU和CPU的计算结果保证了一致。
共享内存的动态分配
我们在数组求和中用到了共享内存
__shared__ float sData[16];
但如果我们定义一个动态数组就会报错，如
int n = 10;
__shared__ float sData[n];
报错如下
main.cu(11): error: expression must have a constant value
main.cu(11): note #2689-D: the value of variable "n"
正确的写法为
#include <iostream>
using namespace std;
__global__ void sum(float *a, float *result) {
int tid = threadIdx.x;
//定义一块同一个线程块的动态共享内存
extern __shared__ float sd[];
float *sData = (float *)sd;
//将每一个线程的值都放入到共享内存中
sData[tid] = a[tid];
//同步所有的线程全部写入共享内存
__syncthreads();
for(int i = 8; i > 0; i /= 2) {
if (tid < i) {
sData[tid] += sData[tid + i];
}
__syncthreads();
}
if (tid == 0) {
result[0] = sData[0];
}
}
int main() {
float a[16];
float result[1];
for(int i = 0; i < 16; i++) {
a[i] = i * (i + 1);
}
float *a_g;
cudaMalloc((void**)&a_g, 16 * sizeof(float));
cudaMemcpy(a_g, a, 16 * sizeof(float), cudaMemcpyHostToDevice);
float *result_g;
cudaMalloc((void**)&result_g, 1 * sizeof(float));
sum<<<1,16,16>>>(a_g, result_g);
cudaMemcpy(result, result_g, 1 * sizeof(float), cudaMemcpyDeviceToHost);
cout << result[0] << endl;
return 0;
}
运行结果
1360
这里共享内存的大小是由
sum<<<1,16,16>>>(a_g, result_g);
中<<<a,b,c>>>的c来分配的。
求两个向量的点积
两个向量的点积为各个分量相乘并相加。
#include <iostream>
using namespace std;
#define LENGTH 16
#define BLOCKNUM 2
#define THREADNUM 4
__global__ void dot(float *a, float *b, float *result) {
int tid = threadIdx.x;
int bid = blockIdx.x;
int total_thread_num = THREADNUM * BLOCKNUM;
__shared__ float sData[THREADNUM];
int global_id = tid + bid * THREADNUM;
sData[tid] = 0;
//将16分成了两段，先算前8个再算后8个
while(global_id < LENGTH) {
sData[tid] += a[global_id] * b[global_id];
global_id += total_thread_num;
}
__syncthreads();
for(int i = THREADNUM / 2; i > 0; i /= 2) {
if (tid < i) {
//计算一个block相加的结果
sData[tid] = sData[tid] + sData[tid + i];
}
__syncthreads();
}
if (tid == 0) {
//获取每一个block的结果
result[bid] = sData[0];
}
}
int main() {
float a[LENGTH];
float b[LENGTH];
float result[BLOCKNUM];
for(int i = 0; i < LENGTH; i++) {
a[i] = i * (i + 1);
b[i] = i * (i - 2);
}
float *a_g;
float *b_g;
float *result_g;
cudaMalloc((void**)&a_g, LENGTH * sizeof(float));
cudaMalloc((void**)&b_g, LENGTH * sizeof(float));
cudaMalloc((void**)&result_g, BLOCKNUM * sizeof(float));
cudaMemcpy(a_g, a, LENGTH * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(b_g, b, LENGTH * sizeof(float), cudaMemcpyHostToDevice);
dot<<<BLOCKNUM, LENGTH>>>(a_g, b_g, result_g);
cudaMemcpy(result, result_g, BLOCKNUM * sizeof(float), cudaMemcpyDeviceToHost);
float res = 0;
//将每一个block的结果最终相加
for(int i = 0; i < BLOCKNUM; i++) {
res += result[i];
}
cout << res << endl;
return 0;
}
运行结果
161432
虽然两个长度16的向量的点积只需要一个block就够了，但我们这里使用了两个block。这里需要注意的是#define定义的是常量而不是变量，所以它可以用在共享内存的分配中。
我们可以用Python来验证一下结果
import numpy as np if __name__ == '__main__': a = np.array([i * (i + 1) for i in range(16)], dtype=np.float32) b = np.array([i * (i - 2) for i in range(16)], dtype=np.float32) c = np.dot(a, b) print(c)
运行结果
161432.0
多维矩阵的内存对齐
以前我们在GPU上分配变量数组的内存空间都是一维的，譬如
cudaMalloc((void**)&a_g, LENGTH * sizeof(float));
但如果是一个二维数组该如何分配呢？
#include <iostream>
using namespace std;
int main() {
size_t width = 10;
size_t height = 10;
//GPU上的二维数组
float *b_g;
//每行宽度可分配的空间大小
size_t pitch;
//分配一段宽度为10，高度也为10的二维数组空间
cudaMallocPitch((void**)&b_g, &pitch, width * sizeof(float), height);
cout << "pitch:" << pitch << endl;
cout << "occupy_num:" << pitch / sizeof(float) << endl;
//释放分配的GPU内存空间
cudaFree(b_g);
return 0;
}
运行结果
pitch:512
occupy_num:128
从结果我们可以看出对于每个二维数组的宽度空间为512，可以存放float类型128个，也就是说我们定义的二维数组的宽度小于128，那么pitch都可以分配512，这里我们可以来试一下
#include <iostream>
using namespace std;
int main() {
size_t width = 128;
size_t height = 10;
float *b_g;
size_t pitch;
cudaMallocPitch((void**)&b_g, &pitch, width * sizeof(float), height);
cout << "pitch:" << pitch << endl;
cout << "occupy_num:" << pitch / sizeof(float) << endl;
cudaFree(b_g);
return 0;
}
打印结果依然是
pitch:512
occupy_num:128
但如果我们将宽度改成129来看一下
#include <iostream>
using namespace std;
int main() {
size_t width = 129;
size_t height = 10;
float *b_g;
size_t pitch;
cudaMallocPitch((void**)&b_g, &pitch, width * sizeof(float), height);
cout << "pitch:" << pitch << endl;
cout << "occupy_num:" << pitch / sizeof(float) << endl;
cudaFree(b_g);
return 0;
}
运行结果
pitch:1024
occupy_num:256
我们可以看到它的宽度空间翻倍了，开启了第二个内存块，变成了1024，可容纳的float也变成了256。它这么做的目的是为了保证二维数组的每行开头地址都是对齐的。
如果是三维数组
#include <iostream>
using namespace std;
int main() {
size_t width = 129;
size_t height = 10;
float *b_g;
size_t pitch;
cudaMallocPitch((void**)&b_g, &pitch, width * sizeof(float), height);
cout << "pitch:" << pitch << endl;
cout << "occupy_num:" << pitch / sizeof(float) << endl;
cudaFree(b_g);
//GPU多维数组的结构体
cudaPitchedPtr c_g;
int x = 10, y = 22, z = 16;
//三维空间大小
cudaExtent extent = make_cudaExtent(sizeof(float) * x, y, z);
//GPU上分配地址空间
cudaMalloc3D(&c_g, extent);
//待拷贝到GPU上的数组
float c[160000] = {0};
//拷贝结构体
cudaMemcpy3DParms cpyParm;
cpyParm = {0};
//待拷贝的数组地址和空间大小
cpyParm.srcPtr = make_cudaPitchedPtr((void*)c, sizeof(float) * width, width, height);
//GPU上的数组地址
cpyParm.dstPtr = c_g;
cpyParm.extent = extent;
cpyParm.kind = cudaMemcpyHostToDevice;
cudaMemcpy3D(&cpyParm);
cudaFree(c_g.ptr);
return 0;
}
这么做并不是唯一的方法，但是它可以帮助我们对齐每行数据的首地址，方便我们更好的操作。
