目录
前言
个人学习目的
主要是用于自己cuda数据挖掘的一些研究和学习，只对重要的部分和数据挖掘相关的内容进行初步的学习，这个文章真的是从0入门
参考资料
使用的书是cuda c编程权威指南，代码pdf网上都有，还有一个视频课https://www.bilibili.com/video/BV11K411N7Wm，视频还附带github，有ppt和代码，这个视频课强烈推荐入门，这里的图片大多引用了这个视频附GitHub的ppt
入门
如何配置环境
首先按https://blog.csdn.net/jhp1266987/article/details/112792853安装好环境，这里是ubuntu系统，使用linux会比较方便，编程语言为cuda，.cuda格式文件，风格是c语言风格
如何运行一个程序
使用nvcc 一些参数 程序.cu -o 输出文件名，执行后再执行./输出文件即可，例如 nvcc arch=sm_60 test.cu -o test，然后执行./test得到输出结果
如何调试
nvprof或者nvvp(nsight），nvprof会输出各部分的运行时间，每秒读写量(GB)等一些信息、nsight会有可视化，具体还没有研究
一些基础知识点
网格grid、块block、warp、thread
大概会生成最高三维度的grid，grid含有多个block，block也最多三维，block包含多个thread，warp是32个thread，是一个虚拟的概念，在规约的时候有用，就不用同步了
具体维度怎么分配影响性能还未知，不过可以通过不断调参调整网格和块的大小，例如2维度的网格，1维度的块；1维度的网格，1维度的块；调整块中线程的个数。不一定哪个性能会好，这个问题还需继续学习原因
SM
SM是流多处理器，具体特点看视频教程，这里只列举个人关心的点，SM包含16或32或64..个bank，大概可以理解为地址，可以处理冲突，例如下左图两个线程要一个bank处理，就是二冲突，右边是8冲突，处理好冲突可以对算法时间进行优化，暂时还没有搞懂其中的奥妙（滑稽）
CUDA执行流程
1.分配host内存，并进行数据初始化;
2.分配device内存，并从host将数据拷贝到device上;
3.调用CUDA的核函数在device上完成指定的运算;
4.将device上的运算结果拷贝到host上 (性能)
5.释放device和host上分配的内存。
host是cpu，device是gpu，大概来讲就是cpu上初始化，gpu初始化，把cpu上需要加速的东西扔给gpu，写个gpu程序执行，然后把结果返回给cpu，释放cpugpu内存
// Kernel定义
__global__ void vec_add(double *x, double *y, double *z, int n)
{
int i = get_tid(); // user-defined function
if (i < n) z[i] = x[i] + y[id];
}
int main()
{
int N = 1000000; // 1M
int bs = 256;
int gs = (N + bs - 1) / bs;
// kernel, call GPU
vec_add<<<gs, bs>>>(x, y, z, N);
}代码大概如上，<<<>>>是调用cuda程序，括号中间是网格和块，__global__是核函数，__global__是异步的，cpu不会等待这个程序会继续执行别的，除非使用一个函数，后面会说
CUDA程序风格和ｃ不同之处
ｃ语言如果是很多数相加的话，可以直接写个ｆｏｒ循环，这里就不用，每个线程都会调用核函数的，核函数中的意思是先获取线程号的范围，然后让所有小于ｎ的号的线程激活，执行下边的加操作
/* get thread id: 1D block and 2D grid */
#define get_tid() (blockDim.x * (blockIdx.x + blockIdx.y * gridDim.x) + threadIdx.x)
/* get block id: 2D grid */
#define get_bid() (blockIdx.x + blockIdx.y * gridDim.x)
__global__ void vec_add(double *x, double *y, double *z, int n)
{
int i = get_tid(); // user-defined function
if (i < n) z[i] = x[i] + y[id];
}内存的使用
有L1，L2，寄存器，这三个是不可编程的，应该说是编程比较费劲
全局内存容量大，但是比共享内存慢，需要对齐操作，就是分成32的倍数，不然会影响性能
共享内存，优化使用，比全局内存快，但是容量不大，只有几十KB
常量内存和纹理内存暂时不需要考虑
规约及优化
规约能做什么，能加法，乘法，逻辑运算
第一种
extern在这里暂时意义不明，__share__是共享内存，可以给share加volatile修饰符，就是禁止硬件自优化，因为自己的优化可能比自动优化更好，所以要禁止自动优化防止降低速度
每个线程把一个元素从全局内存加载到共享内存中，tid和i暂时不知道怎么写
__syncthreads是同步，先等待上边的程序执行完毕，再执行后边的程序，这个标志符号很有用
下边的代码是，每个线程，加两个数，也就是最简单的规约
__global__ void reduce0 (int *g_idata, int *g_odata) {
extern __shared__ int sdata[];
// each thread loads one element from global to shared mem
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
sdata[tid] = g_idata[i];
__syncthreads();
// do reduction in shared mem
for (unsigned int s=1; s < blockDim.x; s *= 2) {
if (tid % (2*s) == 0) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
// write result for this block to global mem
if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}第二种
将下图中，上边的代码替换为下边的，就是不活跃的线程在后续的规约中就不使用了，
for (unsigned int s=1; s < blockDim.x; s *= 2) {
if (tid % (2*s) == 0) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
• With strided index and non-divergent branch
// do reduction in shared mem
for (unsigned int s=1; s < blockDim.x; s *= 2) {
int index = 2 * s * tid;
if (index < blockDim.x / s) {
sdata[index] += sdata[index + s];
}
__syncthreads();
}第三种
这里比较直观，就不说了
// do reduction in shared mem
for (unsigned int s = blockDim.x/2; s > 0; s /= 2) {
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}第四种
一次读取两个数先加起来，就是一种最简单的循环展开
// each thread loads two elements from global to shared mem
// end performs the first step of the reduction
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x* blockDim.x * 2 + threadIdx.x;
sdata[tid] = g_idata[i] + g_idata[i + blockDim.x];
__syncthreads();第五种
剩下最后一个warp。32个数据的时候，直接使用warp的特性相加即可
// do reduction in shared mem
for (unsigned int s = blockDim.x/2; s > 32; s /= 2) {
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
if (tid <32)
{
sdata[tid] += sdata[tid + 32];
sdata[tid] += sdata[tid + 16];
sdata[tid] += sdata[tid + 8];
sdata[tid] += sdata[tid + 4];
sdata[tid] += sdata[tid + 2];
sdata[tid] += sdata[tid + 1];
}最后一种全展开
就是假设block的上限是512个，当然现在的2080ti，3090肯定比这个多，不过处理的方式是一样的，然后全部展开处理，写代码的方式如下
if (blockSize >= 512) {
if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads();
}
if (blockSize >= 256) {
if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads();
}
if (blockSize >= 128) {
if (tid < 64) { sdata[tid] += sdata[tid + 64]; } __syncthreads();
}
if (tid < 32) {
if (blockSize >= 64) sdata[tid] += sdata[tid + 32];
if (blockSize >= 32) sdata[tid] += sdata[tid + 16];
if (blockSize >= 16) sdata[tid] += sdata[tid + 8];
if (blockSize >= 8) sdata[tid] += sdata[tid + 4];
if (blockSize >= 4) sdata[tid] += sdata[tid + 2];
if (blockSize >= 2) sdata[tid] += sdata[tid + 1];
}如果block大小不知道，就用下边的代码，这个代码具有可扩展性，不像上边一下子写死了，不过上限应该不是512，应该再加几行，这个以后再说。
switch (threads)
{
case 512:
reduce5<512><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
case 256:
reduce5<256><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
case 128:
reduce5<128><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
case 64:
reduce5< 64><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
case 32:
reduce5< 32><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
case 16:
reduce5< 16><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
case 8:
reduce5< 8><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
case 4:
reduce5< 4><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
case 2:
reduce5< 2><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
case 1:
reduce5< 1><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
}库函数
cusparse处理稀疏矩阵
cublas线性代数
curand随机数
其他的都用不到，比如openacc，傅里叶变换等
举个例子
在下边的代码是先在gpu上随机10000个数，0或者1，把他们加起来，统计1的个数
首先定义，分配空间，使用generate_uniform_kernel核函数随机数，假设50个块，每个块256个维度，这个是随便设定的，只要50*256>10000即可，如果随机更多的数比如1亿那就要分配更多的块，每个块更多的线程，当然线程数不能超过上限，可以用一个函数查看自己gpu每个块的线程上限，百度就有；
然后执行加法操作，假设5个块256个线程，这样5*256小于10000，所以进入函数后会执行
for (int i = id; i < 10000; i += 1280)
{
shared[tid] += a[i];
__syncthreads();
} id是线程唯一索引，tid是当前块的线程索引每个块都有0-255的线程号，这个一开始我都没搞懂，隔了几天才明白，
shared是共享内存，反正是能加速，别的就不管了，共享内存数量和线程数量一致，好像是共享内存也不能太大，否则影响SM的正常运行，咱也不知道，咱也不敢问。
__synthreads是同步，等待一个块种所有线程做完这个for循环的第一步，然后一起进入循环的第二步，说一起也不是一起，也是按一定顺序来的，但是可能是乱序，这个太细节了，先不细扣
这里id一共是5*256个=1280个
然后把间隔1280的加到一个地方
现在得到了256维度的shared共享内存，要把这256个数加起来，这里用到规约、
传统的规约使用循环，就时间会慢点，这里全部展开，
第一步
if (tid < 128) shared[tid] = shared[tid] + shared[tid + 128];//对循环展开256变成128
__syncthreads();
if (tid < 64) shared[tid] = shared[tid] + shared[tid + 64];//对循环展开128->64
__syncthreads();先把256维度降低到128，再把128降低到64
第二步
if (tid < 32) {//利用了warp的特性，一个warp里就不用同步了，可以省略掉同步的时间
shared[tid] += shared[tid + 32];
shared[tid] += shared[tid + 16];
shared[tid] += shared[tid + 8];
shared[tid] += shared[tid + 4];
shared[tid] += shared[tid + 2];
shared[tid] += shared[tid + 1];
}因为32是一个warp所以不用同步，直接全写出来，这样就可以了。这个优化我从一开始的0.0768ms降低到了0.0168ms，当然线程，块的数量不再多加了，因为越多肯定越快，主要是学习优化的方法
还有一个优化是volatile和unrolling，如果学懂了再更新把，现在学的差不多能写程序了，先做个数据挖掘的baseline出来
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include <curand_kernel.h>
#include<iostream>
#include<thrust/reduce.h>
using namespace std;
const int arraySize = 10000;
__global__ void generate_uniform_kernel(int *a)
{
int id = threadIdx.x + blockIdx.x * blockDim.x;
if (id >= 10000)
return;
curandState state;
curand_init(10, id, 0, &state);
a[id] += curand_uniform(&state) * 2;
//if(id< 100)
//printf("threadIdx,x = %d,blockIdx.x = %d,blockDim.x = %d\n", threadIdx.x, blockIdx.x, blockDim.x);
}
__global__ void sum_c1(int *a, int *c)
{
__shared__ int shared[256];
int id = threadIdx.x + blockIdx.x * blockDim.x;
int tid = threadIdx.x;
if (blockIdx.x == 0)
{
shared[tid] = 0;
}
if (id > 1280)
{
return;
}
for (int i = id; i < 10000; i += 1280)//这个变不了了，因为线程不够，所以数要往前加，如果线程够，就不用这一步骤了
{
shared[tid] += a[i];
__syncthreads();
}
//现在shared有32个数要把他们都加到share[0]里
//int i = blockDim.x / 2;//16
//while (i != 1) {//规约
// if (tid < i) {
// shared[tid] += shared[tid + i];
// }
// __syncthreads();
// i /= 2;
//}
if (tid < 128) shared[tid] = shared[tid] + shared[tid + 128];//对循环展开256变成128
__syncthreads();
if (tid < 64) shared[tid] = shared[tid] + shared[tid + 64];//对循环展开128->64
__syncthreads();
if (tid < 32) {//利用了warp的特性，一个warp里就不用同步了，可以省略掉同步的时间
shared[tid] += shared[tid + 32];
shared[tid] += shared[tid + 16];
shared[tid] += shared[tid + 8];
shared[tid] += shared[tid + 4];
shared[tid] += shared[tid + 2];
shared[tid] += shared[tid + 1];
}
if (tid == 0) {
c[blockIdx.x] =shared[0];//如果id是0，就把结果给result，这是一个块内的运算结果
}
}
int main()
{
int c[arraySize] = { 0 };
int a[arraySize] = { 0 };
int *dev_a = 0;
int *dev_b = 0;
int *dev_c = 0;
cudaSetDevice(0);
// Allocate GPU buffers for three vectors (two input, one output) .
cudaMalloc((void**)&dev_c, arraySize * sizeof(int));
cudaMalloc((void**)&dev_a, arraySize * sizeof(int));
cudaMalloc((void**)&dev_b, arraySize * sizeof(int));
generate_uniform_kernel << <50, 256 >> > (dev_a);
printf("///\n");
cudaEvent_t start, stop;
float elapsedTime = 0.0;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, 0);
sum_c1 << <5, 256 >> > (dev_a, dev_c);
cudaEventRecord(stop, 0);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&elapsedTime, start, stop);
cout <<"time =" <<elapsedTime << "ms"<<endl;
cudaEventDestroy(start);
cudaEventDestroy(stop);
cudaMemcpy(c, dev_c, arraySize * sizeof(int), cudaMemcpyDeviceToHost);
cudaMemcpy(a, dev_a, arraySize * sizeof(int), cudaMemcpyDeviceToHost);
int sum = 0;
for (int i = 0; i < 5; i++)
{
sum += c[i];
}
printf("gpu_sum=%d\n", sum);
sum = 0;
for (int i = 0; i < 10000; i++)
{
sum += a[i];
}
printf("cpu_sum=%d\n", sum);
cudaDeviceReset();
cudaFree(dev_c);
cudaFree(dev_a);
cudaFree(dev_b);
return 0;
}
798
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
