▶ 并行通讯方式：
map 映射 全局一到一 全局单元素计算操作
transpose 转置 一到一 单元素位移
gather 收集 多到一 元素搬运不计算
scatter 分散 一到多 元素搬运不计算
stencil 模板 全局多到一 模板计算（例如卷积）
reduce 归约 全局多到一 元素计算成一个值
scan/sort 扫描排序 全局多到多 元素局部或全局调整
▶ 几种扫描方法：
● 线性扫描，O(n) 个 step（完全不并行），O(n) 次加法。适用于只有一个处理器的情形
[ 1, 2, 3, 4, 5, 6, 7, 8] [ 3 ] // 每个数和它左1格的数字相加 [ 6 ] ... [ 1 3 6 10 15 21 28 36] // 结果
● 闭扫描的 Hillis Steele 算法，O(log n) 个 step（理解为该矩形的宽度），O(n log n) 次加法（理解为该矩形的面积）。适用于处理器较多，算法受步数限制的情形，步骤效率（step efficiency）较高，步数少
[ 1, 2, 3, 4, 5, 6, 7, 8] [ 3 5 7 9 11 13 15] // 每个数和它左1格的数字相加，没有左1格的数字原样补齐 [ 6 10 14 18 22 26] // 每个数和它左2格的数字相加，没有左2格的数字原样补齐 [ 15 21 28 36] // 每个数和它左4格的数字相加，没有左4格的数字原样补齐 [ 1 3 6 10 15 21 28 36] // 结果
● 开扫描的 Blelloch 算法，O(log n) 个 step（HS方法的两倍），O(n) 次加法（把 n 个数字加成一个）。适用于处理器较少，算法受工作量限制的情形，工作效率（work efficiency）较高，步数多
[ 1, 2, 3, 4, 5, 6, 7, 8] [ 3 7 11 15] // 第2k个数和它左1格的数相加，第2k-1个数原样补齐 [ 10 26] // 第4k个数和它左2格的数相加，其他数原样补齐 [ 36] // 第8k个数和它左4格的数相加，其他数原样补齐（其实不需要） [ 10 0] // 写出中间的数（表示原数组前半段的和），最后一个数补0（表示原数组第一个数以前的和） [ 0 10] // 交叉计算，l'=r，r'=l+r [ 3 0 11 10] // 写出前后半段中间的数（表示前后半段各前半段的和） [ 0 3 10 21] // 交叉计算 [ 1 0 3 3 5 10 7 21] // 写出各半段中间的数（表示各半段中前半段的和） [ 0 1 3 6 10 15 21 28] // 交叉计算，结果
▶ compact过程：
input: [ 2, 3, 5, 7,11,13,17,19,23,29,31,37] filter: [ 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1] // 由筛选期计算得到，可以并行 address: [ 0, 1, 1, 1, 2, 3, 4, 5, 5, 5, 6, 6] // filter 的开扫描，若 input[i] 被选中，则它应该放到 output 的第 address[i] 位置 output: [ 2, 7,11,13,17,19,37] // 将输入数组中的值依地址数组相应位置上的值进行输出：[output[address[i]]=input[i] if filter[i] for i in range(len(input))];
▶ APOD: analyze, parallelize, optimice, deploy。
● 并行优化的两种体现：
■ 弱缩放：当并行规模增大时，如何解决更大规模的问题
■ 强缩放：当并行规模增大时，如何缩短解决问题的时间
● CUDA 程序优化方法思路：
■ 优化GPU占用率：各 SM、各 block、各 thread 分配时间相近的任务
■ 合并内存访问，减少全局内存的访问
■ 优化同步延迟：减少 __syncthreads(); 的等待时间，尽量使内存带宽饱和。适量减少每个 block 的 thread 数量，增加每个 SM 的 block 数量来改善，但是若 block 过小、过于分散不利于合并全局内存访问，当线程块数量为流处理器数量的2倍时，计算效率最高（经验关系）
■ 最小化线程分支发散（改进算法）。重整 Warp，使得分支在 Warp 之间而不是 Warp 之内
■ 优化循环结构，减少循环次数差距
■ 使用cuda内置函数，有目的的使用单双精度数字
■ 管理线程通讯（适当增加或减少线程间通讯，调整计算效率）
▶ 7 种常见并行程序优化模式
1、数据布局变换（Data layout transformation）
struct foo { float a; float b; } aa[8]; //结构数组，array of atructures, AOS struct foo { float a[8]; float b[8]; } aa; //数组结构，structure of arrays, SOA
2、发散 - 收集变换（Scatter - to - gather transformation）
out[i]=in[i-10]+in[i]+in[i+10] // 分散地址访问 out[i]=in[i-1]+in[i]+in[i+1] // 紧缩地址访问
3、瓦片（Tiling）利用更广高速度的内存形式，如 __shared__
4、私有化（Privatization）将多个线程需要同时用到的同一内存数据分割或另存为多个副本，供不同线程单独使用，避免内存读取冲突和延迟。例如计算直方图
5、进仓（Binning）将输出位置映射到输入数据的较小子集上。例如筛选距离某点最近的 n 个点，先画粗网格进行第一轮筛选，剩下的点都不要
6、压缩（Compaction）创造一个仅包含活动元素的紧凑数据组，防止多余的线程闲置。加速比例不能超过线程束中的线程数量，即 32 倍。例如大矩阵乘法分块
7、正则化（Regularization）负载均衡，设置每个线程需要完成的工作量的上限，超出的部分利用其它的核函数或者CPU来补充完成。例如典例：地图上寻找相邻的给定点
▶ Warp 含有 32 条 thread，优先按照 threadIdx.x 划分，再按照 threadIdx.y，最后按照 threadIdx.z 划分。同一时刻只能执行统一一条指令，分支结构会先执行一部分，挂起后再执行另一部分，总时间变长。
// 典例： switch(threadIdx.x%32){case 0 ~ 31: foo<<<1,1024>>>();} // 减速为1/32，每个线程分别执行一次 switch(threadIdx.x%64){case 0 ~ 63: foo<<<1,1024>>>();} // 减速为1/32，因为每个Warp中只有32个线程，不可能有更改多的分支 switch(threadIdx.y){case 0 ~ 15: foo<<<1,dim3(64,16)>>>();} // 不减速，因为每个Warp中各线程的threadIdx.y是相等的 switch(threadIdx.y){case 0 ~ 15: foo<<<1,dim3(16,16)>>>();} // 减速为1/2，因为每个Warp中各线程threadIdx.y有两个值 switch(threadIdx.x%2){case 0 ~ 31: foo<<<1,1024>>>();} // 减速为1/2，共 2 种取值 switch(threadIdx.x/32){case 0 ~ 31: foo<<<1,1024>>>();} // 不减速，所有线程计算值相等 switch(threadIdx.x/8){case 0 ~ 31: foo<<<1,1024>>>();} // 减速为1/4，共 4 种取值
▶ 利用宏__CUDA_ARCH__生成同时在主机和设备上运行的同一个程序，并具有不同处理方式，宏__CUDA_ARCH__是一个整数，百位表示计算功能集主版本号
1 __host__ __device__ int myFunc(void) 2 { 3 #if defined(__CUDA_ARCH__) 4 // Device code here 5 #else 6 // Host code here 7 #endif 8 }
▶ GPU工作调度机制：将流中工作映射，先按工作种类（核函数引擎、内存拷贝引擎等）分类再按时间先后串行，不同流的同一类型的工作之间仍然曾在阻塞，应该采用广度优先策略，分拆多个任务穿插道不同的流中，以便在一个流占用核函数引擎的时候另一个流占用内存拷贝引擎
● 若同时运行两个指向同一地址的流，则仍会并行运行，但结果未定义
1 cudaMemcpyAsync(&d_array,&h_array,ARRAY_BYTES,cudaMemcpyHostToDevice,s1); 2 foo<<<blocks,threads,s2>>>(d_array);
▶ 图 G = (V，E) 的并行广度优先遍历算法（O(n2)）：
● 开启V个线程，将根节点标记为0，其他标记为-1
● 开启V个线程，每次循环检查相应的顶点是否存在这样一条边，该边的一端已经被标价，另一端没有被标记：若存在，则将没有被标记的端点标记为已标记的端点的值+1，并且报告遍历尚未结束；若不存在，则不做改变，报告该节点遍历已经结束
1 __global__ void bfs(const Edge * edges, Vertex * vertices, int currentDepth, bool *done) 2 { 3 int e = blockDim.x * blockIdx.x + threadIdx.x; 4 int dfirst = vertices[edges[e].first], dsecond = vertices[edge[e].second]; 5 if (dfirst == currentDepth && dsecond == -1) 6 { 7 vertices[vsecond] = dfirst + 1; 8 *done = false; 9 } 10 else if (dsecond == current_depth && dfirst == -1) 11 { 12 vertices[vfirst] = dsecond + 1; 13 *done = false; 14 } 15 else 16 *done = true; 17 return; 18 }
▶ 图 G = (V，E) 的并行广度优先遍历算法（O(n)）：
● 用类似SCR的方式存储一张图，保存两个数组
C：依次保存每个节点的邻居的编号，长度等于边的条数
R：依次保存每个节点的邻居在 C 中的起点位置，长度等于节点个数+1，最后一个位置存放边的条数，方便最后一个节点的计算
D：依次保存每个节点的深度，长度等于节点个数
● 步骤：
■ 对边界中的每个节点，利用R找到其邻居编号在 C 中的起点以及邻居个数，如对于编号为v的节点，其邻居编号在C中起点为 R[v]，邻居个数为 R[v+1] - R[v]
■ 找到边界的所有相邻节点，依次入队
■ 删除队中已经被标记过的节点（根据D中数据），队空说明已经完成了遍历
■ 确认新节点，标记为边界，返回 1
▶ cuBLAS 使用范例（编译时添加 -L cublas）
1 { 2 int N = 1 << 20; 3 cublasInit(); 4 cublasAlloc(N, sizeof(float), (void **)&d_x); 5 cublasAlloc(N, sizeof(float), (void **)&d_y); 6 7 cublasSetVector(N,sizeof(x[0]), x, 1, d_x, 1); 8 cublasSetVector(N,sizeof(y[0]), y, 1, d_y, 1); 9 10 cublasSaxpy(N, 2.0, x,1, y, 1); // 单精度 y += a*x 11 12 cublasSetVector(N,sizeof(y[0]), d_y, 1, y, 1); 13 14 cublasFree(d_x); 15 cublasFree(d_y); 16 cublasShutdown(); 17 }
▶ MCUDA 工具（Linux平台）将 CUDA 代码编译为可以在主机 CPU 上运行的程序
▶ Thrust库，CUDA中类似STL的并行函数库
▶ CudaDMA库，优化全局内存和共享内存交换
▶ Kahan求和算法：人为记录浮点数加法过程中每一步的舍入误差，并在计算最后进行补偿，减小了总体计算误差
1 { 2 float a[N], temp, compensation = 0.0f, sum_old, sum = 0.0f; 3 for(int i = 0; i < N; i++) 4 { 5 sum_old = sum; //记录前i个数的和 6 temp = a[i] + compemsation; //计算补偿以后的新待加数 7 sum += temp; //获得前i+1个数的和 8 compensation = temp + sum_old - sum0; //计算新的补偿 9 } 10 sum += compensation; //剩余补偿 11 }