cuda并行程序设计复习（直方图、卷积、扫描、前缀和）

程序员文章站 2022-07-12 20:57:03

...

第五章线程执行效率与SIMD

warp线程时单指令多数据执行（SIMD）;warp中的线程执行相同的命令。在任何时间
控制发散：当warp中的线程通过不同的控制决策而采取的不同控制路径时，就会发生，采取不同的控制路径的线程会最终为串行执行，当分支或者循环的条件为线程索引时就可能出现,发生在block中（each block is divided into 32-thread warps
产生的影响取决于数据，对于数据量较大的程序影响较小，对于大型数据而言边界检查所带来的控制发散影响微不足道，并且内核有大量的控制流并不意味着就会有大量的控制发散

第七章直方图计算

直方图: 从大型数据集中国提取出显著特征和模式的方法

例如：图像中的对象识别的特征提取信用卡的交易欺诈检测，天体运动关联

基本直方图
分段分区导致内存访问底下，相邻的吸纳晨光不会访问相邻的内存位置，访问不会合并，DRAM带宽利用率低

交错分区，所有的线程处理一个连续的元素部分，他们都移动到下一个部分并重复，内存访问合并
data race

多个线程对同一块数据进行操作，读取顺序的原因会出现读取的数据经过修改，而丢失一部分操作内容，一共四种类型使用原子操作可以解决
cuda原子操作 atom add sub inc dec min max exch CAS

int atomicAdd(int * address, int val); unsigned int unsigned long long int flaot 四种操作类型

基本直方图代码

__global__ void hist(uchar *b,long size,int *histo)
{
  int i = threadIdx.x+blockIdx.x+blockDim.x;
  int stride = blockDim.x*gridDim.x;
  
 while(i<size){
     int  pos = b[i]-'a';
      if(pos>=0 &&pos <=26){
        atomicAdd(&histo[pos/4],1];
      }
      i+=stride;
  }
}

直方图私有化

创建私有的直方图，会增加副本的空间开销，将所有副本汇总至总和的开销，但是在访问私有化直方图与汇总时，出现竞争的串行情况会减少，性能提升至少10倍

直方图私有化代码

__global void hist_p(uchar*b,long siez,uchar *histo){
  __shared__ int private[7];
  if(theadIdx.x<7)
    private[theadIdx.x] = 0;
  __synctheads();
  
  int i = threadIdx.x + blockIdx.x*blockDim.x;
  int stride = blockDim.x*gridDim.x;
  while(i < size){
    int pos = b[i]-'a';
    if(pos>=0 && pos<=26){
    atomicAdd(&private[pos/4],1);
    }
    i+=stride;
  }
  __syncthreads();
  
  if(threadIdx.x < 7)
    atomicAdd(&histo[threadIdx.x],privat[threadIdx.x]);
}

运算是关联可交换的，直方图的大小应该受限于共享内存，如果太大则可以使用部分直方图私有化的方式，通过范围测试来访问全局内存

八模板运算卷积

卷积：一种数组运算，其中的每一个输出元素都是相邻输入元素的加权和，加权的方式通常由卷积核决定，卷积核在运算时不变

一维卷积代码

__global__ void con_1(flaot *N,float *M,float *P,int Mwid,int wid){
  int i = threadIdx.x + blockIdx.x*blockDim.x;
  
  int start = i - Mwid/2;
  float pValue = 0;
  for(int j=0;j<Mwid;j++){
    if(start+j>=0 && start+j<wid)
      pVaLue += N[start+j] * M[j]
  }
  P[i] = pValue;
}

二维卷积代码

__gloabl__ con_2(uchar *in,uchar* m,uchar *out,int mwid,int w,int h){
  int col = theadIdx.x+blockIdx.x*blockDim.x;
  int row = theadIdx.y+blockIdx.y*blockDim.y;
  
  if(col <w &&row <h){
    int pvalue = 0;
    
    int startc = col - mwid/2;
    int startr = row - mwid/2;
    
    for(int i=0;i<mwid;i++)
      for(int j=0;j<miwd;j++){
          int ci = i+startr;
          int cj = j+startc;
          if(ci>-1 && ci<h && cj>-1 &&cj<w){
            pvalue += in[ci*w + cj] * m[i*mwid+j]
          }
      }
      
      out[row *w +col] =(uchar) pvalue;
  }

tile卷积一维不参与输出

__global__ con1_s(flot* N,float* M,float*P){
  __shared__ float Ns[O_TILE_WID+Mwid-1]
  
  int o_idx = O_TILE_WID*blockIdx.x+threadIdx.x;
  int i_idx = o_idx - Mwid/2;
  if(i_idx>=0 && i_idx<Arraywid){
    Ns[threadIdx.x] = N[i_idx];
  }
  __syncthreads();
  float pvalue = 0;
  if(threadIdx.x < O_TILE_WID){
    for(int j=0;j<Mwid;j++){
      pvalue +=M[j]*Ns[threadIdx.x + j]
    }
    P[o_idx] = pvalue;
    __syncthreads();
  }

}

const __restrict__受限的常量存储，自动适配合适的存储

二维con受限内存

__global__ void con_2(float* P.float* N,int height,int wid,
int channels,const float __restrict__ *M){
  __shared__ float Ns[O_TILE_WID + Mwid-1][O_TILE_WID + Mwid-1]
  int ty = threadIdx.y;
  int tx = threadIdx.x;
  
  int o_col = blockIdx.x*O_TILE_WID + tx;
  int o_row = blockIdx.y*O_TILE_WID + ty;
  int i_col = o_col - Mwid/2;
  int i_row = o_row - Mwid /2;
  if(i_col >=0 && i_col<wid &&i_row >=0 && i_row<height){
    Ns[ty][tx] = N[i_row * pitch +i_col]
  }else  Ns[ty][tx]=0.0f;
  __syncthreads();
  
  float pvalue = 0;
  if(ty<O_TILE_WID && tx <O_TILE_WID){
  
  for(int i=0;i<Mwid;i++)
    for(int j=0;j<Mwid;j++){
      pvalue += M[i][j] * Ns[ty+i][tx+j]
    }
   if(o_col<wid & o_row<heigth)  P[o_row * wid+o_col]=pvalue;
    __syncthreads();
  }
  
}

计算效率
- 对于一维卷积而言，原始的操作需要加载O_TILE_WID*MWID个元素的全局内存访问，而使用啦共享内存后则为 O_TILE_WID+MWID-1个访问次数带宽减少量为两者之比忽略出边界的元素
  - O_TILE_WID*MWID / O_TILE_WID+MWID-1
- 对于二维卷积而言，原始操作带宽较少量为
  - (O_TILE_WID) 2 ^2 2 *MWID 2 ^2 2 / (O_TILE_WID+MWID-1) 2 ^2 2

第九章并行归约

划分与汇总：要求数据集处理的元素没有顺序，关联和交换，可以将数据划分为更小的块，让每一个线程处理一个块，通过归约树将所有分块的结果进行总结求和，类似mapreduce
规约计算：求和求最值阶乘

规约算法的复杂度一般为O(N),N个数值有N次归约操作
归约树算法对于N-1个操作只需要log(N)次，类似二分

对于N个输入值操作N-1次，归约log(N)次。总的并行度为（N-1）/ log(N)

速度上与顺序算法相当，但是资源的使用较大

归约求和每线程处理两个元素到share 内存每个块处理两个dim的大小

__shared__ float Ns[blockDim.x*2]
int t = threadIdx.x
int start = blockDim.x*2*blockIdx.x
Ns[t] = input[start+t]
Ns[blockDim.x+t]= input[start+t+blockDim.x]

for(int stride = 1;stride<blockDim.x;stride *=2){
  __syncthreads(); //在下一步之前需要确保每个元素已经加好
  if(t % stride ==0)
    Ns[2*t] += Ns[2*t+stride];
}

改进的归约

for(int s=blockDim.x;s>0;s=/2){
  __syncthreads();
  if(t < s ){
    p[t]+=p[s+t] 
  }
}

第十章前缀和计算

扫描经常用于并行工作分配和资源分配，基数排序，快速排序，求解递归，

高效的串行算法 O(N)

y[0] = x[0]for i to max_i:  y[i] = y[i-1] +x[i]

 普通的并行 y0 = x0; y1 = x0+x1;...

更有效的前缀和，每个线程加相邻的两个位置的数

__global__ void scan(float* X,float*Y,int inputsize){  __shared__ float XY[SECTIONSIZE]  int i = blockIdx.x*blockDim.x +threadIdx.x;  if(i < inputsize)    XY[threadIdx.x] = X[i];  __syncthreads();  for(int s=1;s<=threadIdx.x,s*=2){     __syncthreads();    float temp = XY[threadIdx.x-s];    __syncthreads();    XY[threadIdx.x] +=temp;  }   __syncthreads();  if(i<inputsize) Y[i] =XY[threadIdx.x];}

工作效率

Total Adds ： n*log(n) -(n-1) => *O(nlog(n))

基于平衡树的前缀和

__global__ void presum(float* go,float* gi,int n){    extern __shared__ float tmp[];    int tid = threadIdx.x;    int offset = 1;        temp[2*tid] = gi[2*tid];    temp[2*tid+1] = gi[2*tid+1]        for(int i = n/2;i>0;i/=2){      __synctheads();      if(tid<i){      int a = offset*(2*tid +1)-1;      int b = offset*(2*tid +2)-1      temp[b] +=temp[a]       }      offset *=2;    }    if(tid==0) temp[n-1] = 0;        for(int i = 1;i<n;i*=2){      offset /=2;      __synctheads();      if(tid < i){          int a = offset*(tid*2+1)-1;          int b = offset*(tid*2+2)-1;                    float t = tmp[a];          tmp[a] = tmp[b];          tmp[b] += t;      }    }  __syncthreads();  go[2*tid] =tmp[tid*2];  go[2*tid+1] = tmp[tid*2+1];}

运行效率

一共**log(n)**次迭代，total Adds n-1 => O(N）

计算归约log(n)-1次迭代 total add n-2 - (log(n)-1) =>O(n)

两者一共不超过 2(N-1) 实际2(n-1)-logn

十一 bank conflict

bank是共享内存上存储单元的划分方式，GPU共享内存是基于这种bank存储体切换的架构
1.x warp=32 bank=16 2.x warp=bank=32 3.x bank可以自定义

每个bank每个周期只能指向一次操作，也就是说每个bank的带宽为每周期 32bit。
什么时候发生？

同一个warp里的不同线程访问同一个bank的不同位置
什么时候不会发生？

广播：当一个warp中的所有线程都访问同一个地址的共享内存时，这是最好的

多播：一个warp中多个线程访问同一个地址的共享内存时，次优

即使同一个warp中的线程随机访问不同的bank，只要没有访问同一个bank的不同位置就不会发生bank conflict
数据并行原则
- 有效利用并行性
- 尽可能合并内存访问
- 使用share memory
- 开发其他内存例如：texture constant
- 减少bank冲突

纹理内存

纹理内存并不是在硬件中对应一块专门的存储器
纹理内存功能：地址映射，数据滤波，缓存等功能
纹理参照系：提供数据与纹理内存的绑定 CUDA Array可以实现对内存，绑定到纹理的线性内存和数组中的元素被称为像元（texels）
texture<Type,Dim,ReadMode>texRef ReadMode读取类型，是否需要归一化默认为否cudaReadModeElementType，是为cudaReadModeNormalizedFloat

struct cudaChannelFormatDesc{ int x,y,z,w;enum cudaChannelFormatKind f;}

cudaChannelFormatKind 指成员类型 cudaChannelFormatKindFloat 浮点型 ~Signed 有符号整型 ~Unsigned 无符号整型

cudaMalloc3DArray() 可以分配1，2，3维的数组，cudaMallocArray() 一般用于二维数据分配，cudaFreeArray() 释放内存，cudaMemcpyToArray()
纹理拾取：纹理拾取函数采用纹理坐标对纹理存储器进行访问。
- 线性内存用texfetch1D函数访问，采用的纹理坐标是整型
- 对与一维、二维和三维CUDA数组绑定的纹理访问: 分别使用tex1D()、tex2D()和tex3D()函数访问，并且使用浮点型纹理坐标。‘

ds_M[ty][tx] = M[G_tx * WIDTH+G_ty + m*TILE]
ds_N[ty][tx] = N[G_ty+m*TILE) * WIDTH+G_tx]

value += ds_M[ty][k] * ds_N[k][tx]
__syncthreads();

P[G_tx*WIDTH+G_ty] = pavlue;

cuda并行程序设计复习（直方图、卷积、扫描、前缀和）

第五章 线程执行效率与SIMD

第七章 直方图计算

八 模板运算 卷积

第九章 并行归约

第十章 前缀和计算

十一 bank conflict

纹理内存

cuda并行程序设计复习（直方图、卷积、扫描、前缀和）

第五章线程执行效率与SIMD

第七章直方图计算

八模板运算卷积

第九章并行归约

第十章前缀和计算