cuda学习

GPU中有多个流处理器SM, 当一个线程块被指定给一个SM后，里面的线程会被划分成线程束（32个线程），在SM上交替运行，也就是说SM上一个时刻只有一个线程束在运行。

函数修饰符：__global__表示该函数只能在GPU上运行，但是可以从CPU或者GPU调用。__device__修饰的函数只能在GPU上运行，并且只能从GPU调用。

cuda 中速度最快是给每个线程私有的寄存器，在核函数中不加修饰声明的变量默认保存在寄存器中，当寄存器不够用时，数据会储存到本地内存中，它实际上和全局内存公用一块储存区域。
加__share__修饰的变量处在在每个块的共享内存中，当此块结束，共享内存被释放。为了避免块内竞争，用__syncthreads()进行同步。

全局内存一般在主机上分配，它的声明周期和应用程序一致。

• 一个线程束由32个连续的线程组成，在一个线程束中，所有的线程按照单指令多线程（SIMT）方式执行；即，所有线程都执行相同的指令，每个线程在私有数据上进行操作。
• 从逻辑角度来看，线程块是线程的集合，它们可以被组织为一维、二维或三维布局。
• 从硬件角度来看，线程块是一维线程束的集合。在线程块中线程被组织成一维布局，每32个连续线程组成一个线程束。
• 一个SM会负责多个ThreadBlock(线程块)的计算任务。每个SP一个时刻负责一个thread。
• 硬件层面，SM中有shared memory, register, L1 cache，因 此ThreadBlock内可以共享shared memory，单独的thread拥有自己的Local memory(先被分配到register中，如果register不够就分配到global memory中)。
• Warp是SM调度和执行的基本单位。SIMT机制使得同一个Warp里的线程根据不同的DATA执行相同的指令。一个SM，一次只能运算一个Block里的一组Warp，如果warp中有线程的DATA没有取到，那么调度下一下warp运算。
• 对齐访问含义就是如果”内存事务“（32和128字节两种）的访问首地址是缓存粒度（L1的128字节或L2缓存的32字节）的偶数倍，即实现了对齐访问。
• 在L1缓存的情况下，由”128字节内存事务“完成访问，如果恰好一个线程束访问的地址是连续的128字节，且首地址恰好又是128的倍数。那么这种访问就说所谓合并访问。
• 共享内存实际上是可受用户控制的一级缓存。申请共享内存后，其内容在每一个用到的block被复制一遍，使得在每个block内，每一个thread都可以访问和操作这块内存，而无法访问其他block内的共享内存。这种机制就使得一个block之内的所有线程可以互相交流和合作。
• 你声明一个二维存储体，编译器会把声明的二维转换成一维线性的，然后再重新整理成二维按照32个存储体，4-Byte/8-Byte宽的内存分布，然后再进行运算的。

规约

template <unsigned int iBlockSize>
__global__ void reduceCompleteUnroll(int * g_idata,int * g_odata,unsigned int n)
{
	//set thread ID
	unsigned int tid = threadIdx.x;
	unsigned int idx = blockDim.x*blockIdx.x*8+threadIdx.x;
	//boundary check
	if (tid >= n) return;
	//convert global data pointer to the
	int *idata = g_idata + blockIdx.x*blockDim.x*8;
	if(idx+7 * blockDim.x<n)
	{
		int a1=g_idata[idx];
		int a2=g_idata[idx+blockDim.x];
		int a3=g_idata[idx+2*blockDim.x];
        int a4=g_idata[idx+3*blockDim.x];
		int a5=g_idata[idx+4*blockDim.x];
		int a6=g_idata[idx+5*blockDim.x];
		int a7=g_idata[idx+6*blockDim.x];
		int a8=g_idata[idx+7*blockDim.x];
		g_idata[idx]=a1+a2+a3+a4+a5+a6+a7+a8;

	}
	__syncthreads();
	//in-place reduction in global memory
	if(iBlockSize>=1024 && tid <512)
		idata[tid]+=idata[tid+512];
	__syncthreads();
	if(iBlockSize>=512 && tid <256)
		idata[tid]+=idata[tid+256];
	__syncthreads();
	if(iBlockSize>=256 && tid <128)
		idata[tid]+=idata[tid+128];
	__syncthreads();
	if(iBlockSize>=128 && tid <64)
		idata[tid]+=idata[tid+64];
	__syncthreads();
	//write result for this block to global mem
	if(tid<32)
	{
		volatile int *vmem = idata;
		vmem[tid]+=vmem[tid+32];
		vmem[tid]+=vmem[tid+16];
		vmem[tid]+=vmem[tid+8];
		vmem[tid]+=vmem[tid+4];
		vmem[tid]+=vmem[tid+2];
		vmem[tid]+=vmem[tid+1];

	}

	if (tid == 0)
		g_odata[blockIdx.x] = idata[0];

}

利用共享内存矩阵转置

__global__ void transformSmemUnrollPad(float * in,float* out,int nx,int ny)
{
    __shared__ float tile[BDIMY*(BDIMX*2+IPAD)];
//1.
    unsigned int ix,iy,transform_in_idx,transform_out_idx;
    ix=threadIdx.x+blockDim.x*blockIdx.x*2;
    iy=threadIdx.y+blockDim.y*blockIdx.y;
    transform_in_idx=iy*nx+ix;
//2.
    unsigned int bidx,irow,icol;
    bidx=threadIdx.y*blockDim.x+threadIdx.x;
    irow=bidx/blockDim.y;  //新矩阵按行读取
    icol=bidx%blockDim.y;
//3.
    unsigned int ix2=blockIdx.y*blockDim.y+icol;
    unsigned int iy2=blockIdx.x*blockDim.x*2+irow;
//4.
    transform_out_idx=iy2*ny+ix2;
    if(ix+blockDim.x<nx&& iy<ny)
    {
        unsigned int row_idx=threadIdx.y*(blockDim.x*2+IPAD)+threadIdx.x;
        tile[row_idx]=in[transform_in_idx];
        tile[row_idx+BDIMX]=in[transform_in_idx+BDIMX];
//5
        __syncthreads();
        unsigned int col_idx=icol*(blockDim.x*2+IPAD)+irow;
        out[transform_out_idx]=tile[col_idx];
        out[transform_out_idx+ny*BDIMX]=tile[col_idx+BDIMX];
 
    }
}

只读内存用于常亮参数读取

__global__ void stencil_1d_readonly(float * in,float * out,const float* __restrict__ dcoef)
{
    __shared__ float smem[BDIM+2*TEMP_RADIO_SIZE];
    int idx=threadIdx.x+blockDim.x*blockIdx.x;
    int sidx=threadIdx.x+TEMP_RADIO_SIZE;
    smem[sidx]=in[idx];

    if (threadIdx.x<TEMP_RADIO_SIZE)

    {
        if(idx>TEMP_RADIO_SIZE)
            smem[sidx-TEMP_RADIO_SIZE]=in[idx-TEMP_RADIO_SIZE];
        if(idx<gridDim.x*blockDim.x-BDIM)
            smem[sidx+BDIM]=in[idx+BDIM];

    }

    __syncthreads();
    if (idx<TEMP_RADIO_SIZE||idx>=gridDim.x*blockDim.x-TEMP_RADIO_SIZE)
        return;
    float temp=.0f;
    #pragma unroll
    for(int i=1;i<=TEMP_RADIO_SIZE;i++)
    {
        temp+=dcoef[i-1]*(smem[sidx+i]-smem[sidx-i]);
    }
    out[idx]=temp;
    //printf("%d:GPU :%lf,\n",idx,temp);
}

cuda锁页内存

对CUDA架构而言，主机端的内存被分为两种，一种是可分页内存（pageable memroy）和页锁定内存（page-lock或 pinned）。可分页内存是由操作系统API malloc()在主机上分配的，页锁定内存是由CUDA函数cudaHostAlloc()在主机内存上分配的，页锁定内存的重要属性是主机的操作系统将不会对这块内存进行分页和交换操作，确保该内存始终驻留在物理内存中。

GPU知道页锁定内存的物理地址，可以通过“直接内存访问（Direct Memory Access，DMA）”技术直接在主机和GPU之间复制数据，速率更快。由于每个页锁定内存都需要分配物理内存，并且这些内存不能交换到磁盘上，所以页锁定内存比使用标准malloc()分配的可分页内存更消耗内存空间

cudaError_t cudaStatus = cudaHostAlloc((void **)&a, size * sizeof(*a), cudaHostAllocDefault);

矩阵相乘

利用pichc操作把不规则矩阵和内存对齐，pitch 长度为新矩阵的列宽(为128倍数)

 float *d_matrix;
 float *dc_matrix = new float[M*N];
 size_t pitch;
cudaMallocPitch(&d_matrix, &pitch, M*sizeof(float), N);
cudaMemcpy2D(d_matrix, pitch, dc_matrix, M* sizeof(float), M * sizeof(float), N, cudaMemcpyHostToDevice);

blockdim.x=256
griddim.x=1000
1000对应a矩阵的行，256对应b矩阵的列(1个线程对应多个列)

__shared__ float a_shared[1000];
for(int j=tix;j<a_cols;j+=bdx){
    a_shared[j]=a[bix*a_pitch+j];
}
__syncthreads();
int i,k;
for(i=tix;i<b_cols;i+=bdx){
    float sum=0;
    for(k=0;k<a_cols;k++){
        sum+=a_shared[k]*b[k*b_pitch+i];
    }
    c[bix*a_pitch+i]=sum;
}

矩阵相乘 (共享矩阵)

__global__ void operator_matmul_h(const float *input1, const float *input2,
                                  float *output, int height, int k, int width,
                                  int broadcast) {
  __shared__ float shared_input1[TILE_SIZE][TILE_SIZE];
  __shared__ float shared_input2[TILE_SIZE][TILE_SIZE];

  int batch_idx = blockIdx.z;
  if (broadcast != 1) input1 += batch_idx * height * k;
  if (broadcast != 2) input2 += batch_idx * k * width;
  output += batch_idx * height * width;

  int bx = blockIdx.y;
  int by = blockIdx.x;
  int tx = threadIdx.y;
  int ty = threadIdx.x;

  int row = bx * TILE_SIZE + tx; // x代表行 ，y代表列，这里反过来了
  int col = by * TILE_SIZE + ty;
  float v = 0;

  for (int i = 0; i < (int)(ceil((float)k / TILE_SIZE)); i++) {
    if (i * TILE_SIZE + ty < k && row < height)
      shared_input1[tx][ty] = input1[row * k + i * TILE_SIZE + ty];
    else
      shared_input1[tx][ty] = 0;

    if (i * TILE_SIZE + tx < k && col < width)
      shared_input2[tx][ty] = input2[(i * TILE_SIZE + tx) * width + col];
    else
      shared_input2[tx][ty] = 0;
    __syncthreads(); // k/TILE_SIZE个子矩阵中的第i个拷贝完毕

    for (int j = 0; j < TILE_SIZE; j++)
      v += shared_input1[tx][j] * shared_input2[j][ty];
    __syncthreads(); // 第 i个子矩阵的[tx][ty]值更新
  }

  if (row < height && col < width) output[row * width + col] = v;
  // k/TILE_SIZE个子矩阵个子矩阵的[tx][ty]位置的和为最终结果中[row][width]的值
}

stream

• CUDA流表示一个GPU操作队列，该队列中的操作将以添加到流中的先后顺序而依次执行。可以将一个流看做是GPU上的一个任务，不同任务可以并行执行。使用CUDA流，首先要选择一个支持设备重叠（Device Overlap）功能的设备，支持设备重叠功能的GPU能够在执行一个CUDA核函数的同时，还能在主机和设备之间执行复制数据操作。
  支持重叠功能的设备的这一特性很重要，可以在一定程度上提升GPU程序的执行效率。一般情况下，CPU内存远大于GPU内存，对于数据量比较大的情况，不可能把CPU缓冲区中的数据一次性传输给GPU，需要分块传输，如果能够在分块传输的同时，GPU也在执行核函数运算，这样的异步操作，就用到设备的重叠功能，能够提高运算性能

• 首先声明一个Stream，可以把不同的操作放到Stream内，按照放入的先后顺序执行。
  cudaMemcpyAsync操作只是一个请求，表示在流中执行一次内存复制操作，并不能确保cudaMemcpyAsync函数返回时已经启动了复制动作，更不能确定复制操作是否已经执行完成，可以确定的是放入流中的这个复制动作一定是在其后 放入流中的其他动作之前完成的。

#include "cuda_runtime.h"  
#include <iostream>
#include <stdio.h>  
#include <math.h>  
 
#define N (1024*1024)  
#define FULL_DATA_SIZE N*20  
 
__global__ void kernel(int* a, int *b, int*c)
{
	int threadID = blockIdx.x * blockDim.x + threadIdx.x;
 
	if (threadID < N)
	{
		c[threadID] = (a[threadID] + b[threadID]) / 2;
	}
}
 
int main()
{
	//获取设备属性
	cudaDeviceProp prop;
	int deviceID;
	cudaGetDevice(&deviceID);
	cudaGetDeviceProperties(&prop, deviceID);
 
	//检查设备是否支持重叠功能
	if (!prop.deviceOverlap)
	{
		printf("No device will handle overlaps. so no speed up from stream.\n");
		return 0;
	}
 
	//启动计时器
	cudaEvent_t start, stop;
	float elapsedTime;
	cudaEventCreate(&start);
	cudaEventCreate(&stop);
	cudaEventRecord(start, 0);
 
	//创建一个CUDA流
	cudaStream_t stream;
	cudaStreamCreate(&stream);
 
	int *host_a, *host_b, *host_c;
	int *dev_a, *dev_b, *dev_c;
 
	//在GPU上分配内存
	cudaMalloc((void**)&dev_a, N * sizeof(int));
	cudaMalloc((void**)&dev_b, N * sizeof(int));
	cudaMalloc((void**)&dev_c, N * sizeof(int));
 
	//在CPU上分配页锁定内存
	cudaHostAlloc((void**)&host_a, FULL_DATA_SIZE * sizeof(int), cudaHostAllocDefault);
	cudaHostAlloc((void**)&host_b, FULL_DATA_SIZE * sizeof(int), cudaHostAllocDefault);
	cudaHostAlloc((void**)&host_c, FULL_DATA_SIZE * sizeof(int), cudaHostAllocDefault);
 
	//主机上的内存赋值
	for (int i = 0; i < FULL_DATA_SIZE; i++)
	{
		host_a[i] = i;
		host_b[i] = FULL_DATA_SIZE - i;
	}
 
	for (int i = 0; i < FULL_DATA_SIZE; i += N)
	{
		cudaMemcpyAsync(dev_a, host_a + i, N * sizeof(int), cudaMemcpyHostToDevice, stream);
		cudaMemcpyAsync(dev_b, host_b + i, N * sizeof(int), cudaMemcpyHostToDevice, stream);
 
		kernel << <N / 1024, 1024, 0, stream >> > (dev_a, dev_b, dev_c);
 
		cudaMemcpyAsync(host_c + i, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost, stream);
	}
 
	// wait until gpu execution finish  
	cudaStreamSynchronize(stream);
 
	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop);
 
	std::cout << "消耗时间： " << elapsedTime << std::endl;
 
	//输出前10个结果
	for (int i = 0; i < 10; i++)
	{
		std::cout << host_c[i] << std::endl;
	}
 
	getchar();
 
	// free stream and mem  
	cudaFreeHost(host_a);
	cudaFreeHost(host_b);
	cudaFreeHost(host_c);
 
	cudaFree(dev_a);
	cudaFree(dev_b);
	cudaFree(dev_c);
 
	cudaStreamDestroy(stream);
	return 0;
}