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

2021-12-29 21:00:04 阅读：172 来源： 互联网

标签：__ 并行程序 float int 直方图线程 cuda TILE 内存

第五章线程执行效率与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;

标签：__,并行程序,float,int,直方图,线程,cuda,TILE,内存
来源： https://blog.csdn.net/weixin_46439063/article/details/122224450

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

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

第五章 线程执行效率与SIMD

第七章 直方图计算

八 模板运算 卷积

第九章 并行归约