CUDA流总结

1. 流常用的基础知识

1. 什么是流：CUDA的所有操作均位于一个流中，可将流视为一条流水线，cuda操作视为流水线上的工作
2. 同步与异步（重点，且容易出bug）
   同步：阻塞主机端线程，直到完成
   异步：被调用后，控制权立即返回主机端，无需等待完成

cudaMemcpy为同步，cudaMemcpyAsync为异步

3. 空流与非空流
   空流：也是默认流，不需要显示声明，由于没有ID不利于管理
   非空流：通过cudaStreamCreate显示声明

空流是同步流，非空流是异步流

4. 阻塞与非阻塞
   阻塞流中的操作会被阻塞
   示例：

		kernal_1<<<grid, block,0, stream[0]>>>();
		kernal_2<<<grid, block>>>();
		kernal_3<<<grid, block,0,stream[1]>>>();

上述代码中，kernal_3被kernal_2阻塞，2等待1执行完才会在GPU上执行，3等2执行完再执行

2. 流的一般使用步骤

1. 声明流
2. 分配空间
3. 为每个流分配任务

示例：

    cudaStream_t stream[n_stream];
    for(int i=0;i<n_stream;i++)
    {
        CHECK(cudaStreamCreate(&stream[i]));
    }
       for(int i=0;i<n_stream;i++)
    {
        int ioffset=i*iElem;
        int blockoffset=i*blocksPerGrid/n_stream;
        cudaMemcpyAsync(&a_d[ioffset],&a_h[ioffset],nByte/n_stream,cudaMemcpyHostToDevice,stream[i]));
        findMaxOfMatrix<<<grid,block,0,stream[i]>>>(&a_d[ioffset],&maxPerBlocks_d[blockoffset]);
        cudaMemcpyAsync(&maxPerBlocks_h[blockoffset],&maxPerBlocks_d[blockoffset],sizeof(int)*(blocksPerGrid/n_stream),cudaMemcpyDeviceToHost,stream[i]);
    }

使用流的关键点

1. 同步与异步的控制：异步操作不会等到执行完毕后再执行下一条指令。因此，如果某条指令需要使用上一条指令的结果，此时如果上一条指令是异步的，可能会出现错误
2. 数据量的切分与地址偏移
   流操作经常将一份大数据切分成很多小数据块，通过并行执行小数据块实现隐藏延迟。因此每个流处理的数据块是不同的，需要通过移动指针来控制。

3. 与OpenMP的联合使用

1. omp_set_num_threads(n_stream);
   指定OpenMP并行区域需要调用的CPU核心数
2. #pragma omp parallel
   将后面的花括号内标记为并行部分
3. int i=omp_get_thread_num();
   建立映射，每个主机线程返回唯一一个线程ID

一般使用模板

	omp_set_num_threads(n_stream);
	#pragma omp parallel
	{
		int i=omp_get_thread_num();
		kernal_1<<<grid, block>>>();
		kernal_2<<<grid, block>>>();
		...
	}

4. 流回调

本质是一个函数，可在流回调之前所有流操作执行完毕后调用
使用方法：

//先定义一个自己的流回调函数
void CUDART_CB my_callback(cudaStream_t stream, cudaError_t status, void *data){ ...}

//在流中期望的位置添加下列函数
cudaStreamAddCallback(stream[i], my_callback, (void *)(stream_ids+i), 0);

由于笔者在优化寻找最大值案例时，使用流回调反而导致了耗时增加，因此后面案例中不涉及流回调的使用。

5. 案例：使用流在给定数组内寻找最大值的功能

/*
	Author：Zheng Yeping
	Data:2021-7-10
	Description: using streams to find max
*/

#include <cuda_runtime.h>
#include <stdio.h>
#include <time.h>
#include <windows.h>
#define CHECK(call)\
{\
  const cudaError_t error=call;\
  if(error!=cudaSuccess)\
  {\
      printf("ERROR: %s:%d,",__FILE__,__LINE__);\
      printf("code:%d,reason:%s\n",error,cudaGetErrorString(error));\
      exit(1);\
  }\
}



#define Height  512
#define Width   1024 
#define n_stream 4
const int N = Height * Width;
const int threadsPerBlock = 512;
const int blocksPerGrid = 1024;

void initDevice(int devNum)
{
  int dev = devNum;
  cudaDeviceProp deviceProp;
  CHECK(cudaGetDeviceProperties(&deviceProp,dev));
  printf("Using device %d: %s\n",dev,deviceProp.name);
  CHECK(cudaSetDevice(dev));
}

__global__ void findMaxOfMatrix(int* d_a, int* maxPerBlocks)
{
	__shared__ int partialMax[threadsPerBlock];
	int i = threadIdx.x + blockIdx.x * blockDim.x;
	int tid = threadIdx.x;

	partialMax[tid] = d_a[i / Width * Height + i % Width];

	__syncthreads();
	if (blockDim.x>=512 && tid < 256)
	{
		if (partialMax[tid] < partialMax[tid + 256])
			partialMax[tid] = partialMax[tid + 256];

	}
	__syncthreads();
		if (blockDim.x>=256 && tid < 128)
	{
		if (partialMax[tid] < partialMax[tid + 128])
			partialMax[tid] = partialMax[tid + 128];

	}
	__syncthreads();
		if (blockDim.x>=128 && tid < 64)
	{
		if (partialMax[tid] < partialMax[tid + 64])
			partialMax[tid] = partialMax[tid + 64];

	}
	__syncthreads();

	if (tid < 32)
	{
		if (partialMax[tid] < partialMax[tid + 32])
			partialMax[tid] = partialMax[tid + 32];
		if (partialMax[tid] < partialMax[tid + 16])
			partialMax[tid] = partialMax[tid + 16];
		if (partialMax[tid] < partialMax[tid + 8])
			partialMax[tid] = partialMax[tid + 8];
		if (partialMax[tid] < partialMax[tid + 4])
			partialMax[tid] = partialMax[tid + 4];
		if (partialMax[tid] < partialMax[tid + 2])
			partialMax[tid] = partialMax[tid + 2];
		if (partialMax[tid] < partialMax[tid + 1])
			partialMax[tid] = partialMax[tid + 1];
	}
__syncthreads();
	if (tid == 0)
		maxPerBlocks[blockIdx.x] = partialMax[0];
}

int main(int argc,char **argv)
{
    // set up device
    initDevice(0);
    float time1;
    int nElem=N;
    int nByte=sizeof(int)*nElem;
    int * a_h,*maxPerBlocks_h;
    CHECK(cudaHostAlloc((int**)&a_h,nByte,cudaHostAllocDefault));
    CHECK(cudaHostAlloc((int**)&maxPerBlocks_h,blocksPerGrid*sizeof(int),cudaHostAllocDefault));
    
   cudaMemset(a_h,0,nByte);
   cudaMemset(maxPerBlocks_h,0,blocksPerGrid*sizeof(int));

    int *a_d,*maxPerBlocks_d;
    CHECK(cudaMalloc((int**)&a_d,nByte));
    CHECK(cudaMalloc((int**)&maxPerBlocks_d,blocksPerGrid*sizeof(int)));

    int maxval_cpu=0;
	//generate a random Mat
	for (int i = 0; i < Height*Width;i++)
		{
			a_h[i] = (rand() % 2018) + 1;
			if (maxval_cpu<a_h[i])
			{
				maxval_cpu = a_h[i];
			}
		}

    dim3 block(threadsPerBlock);
    dim3 grid(blocksPerGrid/n_stream);
    
    //asynchronous calculation
    int iElem=nElem/n_stream;
    cudaStream_t stream[n_stream];
    for(int i=0;i<n_stream;i++)
    {
        CHECK(cudaStreamCreate(&stream[i]));
    }
    cudaEvent_t start,stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start,0);
    for(int i=0;i<n_stream;i++)
    {
        int ioffset=i*iElem;
        int blockoffset=i*blocksPerGrid/n_stream;
        CHECK(cudaMemcpyAsync(&a_d[ioffset],&a_h[ioffset],nByte/n_stream,cudaMemcpyHostToDevice,stream[i]));
        findMaxOfMatrix<<<grid,block,0,stream[i]>>>(&a_d[ioffset],&maxPerBlocks_d[blockoffset]);
        CHECK(cudaMemcpyAsync(&maxPerBlocks_h[blockoffset],&maxPerBlocks_d[blockoffset],sizeof(int)*(blocksPerGrid/n_stream),cudaMemcpyDeviceToHost,stream[i]));
    }

    cudaMemcpy(&maxPerBlocks_h, &maxPerBlocks_d, sizeof(int)*blocksPerGrid,cudaMemcpyDeviceToHost);
	//cudaDeviceSynchronize();
    //改用cudaMemcpy进行隐式同步，比显式同步更省时间

    CHECK(cudaEventRecord(stop, 0));
    CHECK(cudaEventSynchronize(stop));

	//如果将下面一段代码放到事件函数上面，结果会出错 :)
	int max = 0;
	for (int i = 0; i < blocksPerGrid; i++)
		if (max < maxPerBlocks_h[i])
			max = maxPerBlocks_h[i];

  	cudaEventElapsedTime(&time1, start, stop);
    printf("\n The time of  calculating is : %f ms\n", time1);
 	printf("\n CPU_maxNum:%5d ", maxval_cpu);
	printf("\n GPU_maxNum:%5d ", max);
    for(int i=0;i<n_stream;i++)
    {
        CHECK(cudaStreamDestroy(stream[i]));
    }
    cudaFree(a_d);
    cudaFree(maxPerBlocks_d);
    cudaFree(a_h);
    cudaFree(maxPerBlocks_h);
    return 0;
}