拒绝访问 | www.chiphell.com | 百度云加速
请打开cookies.
此网站 (www.chiphell.com) 的管理员禁止了您的访问。原因是您的访问包含了非浏览器特征(431aafb-ua98).
重新安装浏览器，或使用别的浏览器随便写写，随便看看
CUDA加速计算项目采坑总结
CUDA加速运算项目采坑总结
本篇博客主要记录研一上期末考试后，完成CUDA加速计算项目(C++)的采坑旅程。
1.Nsight5.2 BUG 不支持gtx1050和gtx1050Ti
Nsight是N厂开发用来配合Visual Studio调试GPU（device端）内存的插件。这个BUG，NVIDIA负全责。Nsight一般会随cuda的安装配套安装（我的cuda是8.0，配套Nsight5.2），但是坑爹的地方来了。cuda8.0有两个版本（具体可以在NVIDIA官网查到），所以Nsight5.2也是两个版本，官网的User Guide显示5.2是支持gtx1050和gtx1050Ti的，但是我的5.2是不支持的，只支持上半年的，1080。
解决方法：
升级Nsight版本。
2.GPU浮点数计算精度问题
GPU端计算浮点数精度和CPU端不同。GPU默认是float精度，CPU默认是double精度。这就意味着在GPU作浮点数乘除计算结果如果累计的话，必然存在误差，并且误差会随累计从10^-6逐渐变大。
解决方法：
强制类型转换：ans =(double)a*(double)b;消除浮点数计算的累计误差。
3.kernel内部申请动态内存问题
官方文档有说明，计算能力在2.0以上的显卡都可以在kernel内部申请动态内存。语法同CPU端C/C++，malloc和new都可以。但是一定要记得free和delete。不过这样做分配内存的时间效率比在CPU（host端）用cudaMalloc之类的函数效率要低。另外就是动态内存大小也有上限，因为thread多了之后还蛮容易溢出的，需要注意。
这部分具体的可以参考。
解决方法：动态内存申请了记得释放。
4.kernel内部计算没有被执行，感觉像被直接跳过了
这是因为块内的寄存器（register）不够了。你的显卡一个块内具体可以有多少寄存器可以通过launch bonds去计算，好像也有相关命令去改变cuda默认限制。
解决方法：
其实首先可以考虑改变整个程序架构，是不是设计的有问题，因为一般是不会register不足的。不然就该大限制或者优化变量内存。
5.unspecified launch failure in XXX 具体代码函数
这个BUG就有意思了，通俗来讲就是说GPU他平常是CPU的小弟，给老板打工的，CPU每隔一段时间就要压榨一下GPU让他来更新一下屏幕，具体时间间隔看系统（Win10是2s，xp和linux都是5s）。如果你现在压榨GPU给你打工（CUDA做计算）超过这个时间，CPU老板就会不高兴，就会报错。这个机制专业术语叫TimeoutDetection and Recovery mechanism，TDR。
解决方法：
Windows系统下可以参考，这里感谢知乎用户：Russell 提供帮助。
没有更多推荐了，
加入CSDN，享受更精准的内容推荐，与500万程序员共同成长！为什么需要GPU
几年前我启动并主导了一个项目，当时还在谷歌，这个项目叫谷歌大脑。该项目利用谷歌的计算基础设施来构建神经网络。
规模大概比之前的神经网络扩大了一百倍，我们的方法是用约一千台电脑。这确实使深度学习取得了相当大的进展。用到相当多的
计算机。不久之后我发现，之前我并没意识到，用一千台电脑是一项非常昂贵的技术。因此，我和我的朋友，意识到，利用一种
不同的技术，仅用三台电脑，而非一千台，就可以做到这点，而秘诀就是利用GPU技术。&
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ---Andrew Ng
GPGPU是最近几年才流行起来的Project，并且在科学计算上使用频繁。
深度学习在06年立项，当时的数据集小的可怜，只需要CPU仿真一下。
09年，Stanford的CV实验室公布ImageNet后，深度学习开始拥抱大数据。
浅层结构+大数据尚且依赖云计算，深度结构+大数据的需求则是云计算提供不了的。
12年，Hinton组的Alex开源了cuda_convnet，正式大规模使用GPU。衍生出Theano、Caffe等项目。
在今天，如果你要修改这些项目代码，包括自己的创新改良，不会CUDA是不行的。
这也是深度学习带来的一个全新领域，它要求研究者不仅要理论强，建模强，程序设计能力也要过硬，不能纸上谈兵。
单从这点来看，深度学习依然和传统机器学习学派分道而行（多使用CPU多线程、云计算）。
大规模分布式CPU云应该尽量避免，尽管面子大，唬人效果好，但毕竟计算效率低，平台组建麻烦。
高效的程序设计，应当从单机中，榨取更多的计算力，进而考虑多机联合。
CUDA物理计算体系
1.1 CUDA特性
CUDA可以说是介于NVIDIA Driver和DirectX等库的中间产物。
传统游戏开发者，无须考虑太多GPU硬件架构，毕竟A卡、N卡水火不容。
通过设备商提供的Driver，DirectX或者OpenCL可以将不同架构的硬件底层封装，提供给开发者共有的API。
而CUDA，则可以看成是Driver、DX之流的上级。
为了泛型编程(C、C++、Fortran多语言）、以及榨取更多的计算力，NVIDIA对OpenCL进行的改装，
贴合自己的GPU硬件架构，量身定做出CUDA。
较游戏程序员不同，CUDA程序员主要工作，就是把握硬件架构，在算法理论时间复杂度下，
将算法串行执行体系，改组为并行执行体系。违背硬件的并行设计，只会让你的程序跑的百倍之慢。
1.2 NVIDIA GPU三大经典架构（Fermi、Kepler、Maxwell)
与海贼王里CP9 (道力可以反映一个人的战斗力，但不能完全决定一个人的实力)设定相似：
NVIDIA为它在不同成长阶段卖出的产品，规定了计算能力体系：
计算能力1.0是跑CUDA的最低条件，这一时期的代表作是8800GT家族。
Fermi架构的计算能力是2.0, Kepler是3.0，Maxwell是4.0。
1.2.1 年轻有为的Fermi（费米）
Fermi架构，也就是GTX 4XX系卡，这算是CUDA的一个重要转折点&&它把流处理器（SP）改名为"CUDA核心"。
NVIDIA的：
NVIDIA的物理计算体系以分为SM、CUDA核心两大部分。
官方称SM为Stream Multiprocessors（流多处理器），民间多译为SM计算阵列。
物理计算体系对于程序员是透明的，需要调度NVIDIA提供的逻辑计算体系（线程Grid、线程Block、线程Thread)
对于一个指定好逻辑计算体系的任务，GPU会按照负载均衡的原则，将任务平均分至各个SM阵列，
对于每个SM阵列，就调度它手下那一伙CUDA核心干活。
由于每个SM阵列的CUDA核心有限，NVIDIA规定，Fermi架构，每个SM最多并行执行1024个线程。
当然，实际任务中，每个SM会分到几百万个线程，这时候，就只能小部分并行，然后再串行了。
Fermi 1.0架构，官方设计是16组SM，512SP，然而旗舰GTX480最后只弄出了15组，480SP，顺次阉割出GTX470、GTX460。
Fermi 2.0架构，旗舰GTX580，总算达到设计图要求，达到了16组，512SP，顺次阉割出了GTX570，GTX560。
1.2.2 老当益壮的Kepler（开普勒）
Kepler架构最大变化在于， 对每个SM阵列，将SP数量扩大到6倍，达到192SP。谓之曰SMX阵列。
每个SMX阵列，包含192个CUDA核心，单次并行吞吐量是2048个线程。
Kepler&1.0架构，官方设计是15组SM，2880SP，然而旗舰GTX580最后只弄出了8组，1536SP，顺次阉割出GTX570、GTX560。
Kepler 2.0架构，旗舰GTX680，总算达到设计图要求，达到了15组，2880SP，顺次阉割出了GTX670，GTX660。
特别版，GTX Titan Z，直接把两块GK110并在一起，合出了30组，5760SP，同时支持双精度浮点计算。
其阉割掉双精度之后，就是GTX690。
值得一提的是，GTX游戏卡直接把双精度阉割掉了，因为只有Tesla做科学计算的时候，才会用双精度浮点运算。
整体来看，Kepler并没有提升SM阵列数量（在产量情况下，GTX 5XX居然只造出了8组）。但是每个阵列的吞吐量扩大到2倍（）。
速度也较之提升（CUDA核心扩大到6倍）。
在GTX 6XX系中，15组SM阵列，总算达到了Fermi的两倍吞吐量。
1.2.3 "长者"风范的Maxwell（麦克斯韦）
Maxwell架构是老黄的无奈之举。因为台积电把20nm工艺让给了ARM系（Apple和高通）。
还是基于28nm的Maxwell，继续在SM上大刀阔斧闹改革，将192SP降低为128SP，谓之曰SMM阵列。
Maxwell 最新架构，官方设计是16组SM，2048SP，为旗舰GTX980，顺次阉割出GTX970，GTX960。
特别版，GTX TitanX，24组SM，3072SP，较之TitanZ，阉掉了双精度浮点数支持。
TitanX是老黄在GTC 2015向DL界主推的一块民用卡，因为DL无需高精度浮点，用Tesla太奢侈。
(仔细分析老黄那自信满满进攻DL领域的表情，拿TitanX打游戏实在有点.....)
降低了SP数量之后，Maxwell得到了功耗和性能之间的权衡调整。
1.3 设备信息
编译CUDA Examples下的deviceQuery，你可以看到GPU的基本信息。
红字标出的，是设计并行程序时需要参考的几个重要设备参数。
CUDA逻辑计算体系
2.1 线程网格(Grid)、线程块（Block)、线程(Thread)、线程束（Warp)
2.1.1 内核函数
内核函数是并行计算中最基本的单元函数，其特点是：
统一的处理逻辑代码，分布并行掌控不同区域的数据，以此达到多区域数据联动并行执行。
NVIDIA为了CPU在逻辑上，能调度GPU计算的函数，规定了统一的格式。
以__global__限定符为始，声明：__global__ void helloworld()。
__global__意思为，GPU执行，CPU调用
调用时，需要分配& &&&线程块，块内线程数&&&。
如执行helloword，使用1个线程块，块内使用256个线程，则
helloworld&&&1,256&&&
2.1.2 线程网格(Grid)
线程网格在编程时并不存在，它只是抽象上的并行网格体系。
不同种类的内核函数，每种内核函数调度数个的线程块，这数个线程块逻辑上被判为一个Grid。
2.1.3 线程块(Block)
线程块是一个3D结构，强调3D坐标系时，需要以dim3类型声明三维大小。
dim3是个结构体, 成员x、y、z，代表方向轴尺度。
如helloworld&&&dim3(1,1,1),256&&&。
当然，大部分操作基本使用的是1D坐标系，线程块默认全部扩展到X轴上。
一般写成helloworld&&&1,256&&&。
通常在内核函数内，需要获取线程块编号，以便对数据集的不同区域处理，四大重要属性：
?dim3 gridDim(不是指有多少Grid，而是指一个Grid有多少Block)
?dim3 blockDim(不是指有多少Block，而是指一个Block有多少Thread)
?dim3 blockIdx
?dim3 threadIdx
对于1D坐标系，有int tid= (blockDim.x*blockIdx.x) + threadIdx.x;
tid指明当前线程的编号，是内核函数里最基本的控制变量。
int step=(blockDim.x * gridDim.x);
由于CUDA限制每个Block的线程数(2.0以上通常使用1024，以下通常使用512）
所以在常规元素分解模型中，通常把每个Block的线程数设置为常量(固定不动)
这时，有两个策略：
① 其一，不固定Block：
这种方法最为常用，由于CUDA对每个任务而言，对Block数量的限制很松，
这时候，可以采取为每个线程分配一个元素的方法，用
$BLOCKS=(N + THREADS - 1) / THREADS$
算出一个动态的Block数量的需求，这时候，for(i=i&N;i+=step)等效于for(i=i&N;i+=1)
因为这个循环根本不会执行第二次。
① 其二，固定Block数量：
这时，这时候，为了跑完全部的N个元素，有些Thread会启动人工循环。
i+=step会将元素坐标继续跳转，因为N必然大于step，你不能用+1来取剩余的元素吧？
这两种方法本质上是等效的，由于在物理执行时，同时并行线程最多大概是3072，
几百万、甚至几千万的Block会被CUDA扔到等待队列里，由CUDA自己安排自动循环，
有时甚至比你的人工循环更高效，所以，通常用①的方法，为了保持写法一致，step也会作为默认跳转量。
2.1.4 线程(Thread)
CUDA逻辑体系里最基本的执行单位，等效于CPU的线程。
内核函数一旦被指明了线程块大小，线程大小后，每个线程就分配到了一个内核函数的副本。
区别这些的线程的唯一方法就是线程编号tid，通过tid，让不同线程窥视数据集的不同部分。
用相同的逻辑代码，执行数据空间的不同子集。
2.1.5 线程束(Warp)
线程束对用户透明，它是NVIDIA强行规定的。目前显卡都固定为32。
逻辑上，线程束将32个线程编为一组。
一般微机系统，如8086，它的访存模式是串行的。每一个总线周期，吞一个字节进来。
NVIDIA的GPU在一个总线周期内，能够最大吞32*4=128字节。
前提是当个线程束内的线程，逐序访问显存，这特别需要设计数据存储形式。
使用线程束的目的是掩盖单个总线周期过长的问题，通常要跑500~600个T周期。
简单样例: HelloWorld、向量加法
3.1 CUDA程序基本结构
CUDA的编译器NVCC只能编译.cu文件。cu文件可以看作是h&cpp文件的混合体。
所有代码可以写在一个cu文件里，但是毫无可读性。
CUDA最大的亮点在于提供C(C++)接口，所以设计模式应该保持传统C程序风格。
?风格一：将不同函数写进不同cu文件里，模块化。
看起来挺好的，实际操作起来比较麻烦。
NVCC编译器在生成程序时，不支持跨文件link，也就是说，#include一个cu文件，会导致重定义。
VS中会这么报错：
error LNK2005: "void __cdecl cudaStrcpy(char *,char *)" (?cudaStrcpy@@YAXPAD0@Z) 已经在 gpu_helloworld.cu.obj 中定义
所以，NVCC编译完cu文件后，要在项目属性里把冲突的cu文件从生成中排除，然后link。
?风格二：h、cpp、cu混用
上图中指明了，h&cpp是由C编译器来编译的。C编译器里不允许#include一个cu文件（不支持)
否则VS会报错
main.obj : error LNK2001: 无法解析的外部符号 _threadIdx
解决方案：
若要引用cu里的函数，在main.cpp里外部extern声明一下，让VS转为NVCC编译器处理。
3.2 GPU版HelloWorld
3.2.1 内存与显存之争
&&老板，我要4G显存的显卡。
&&& &&亲，这是我们最新款显卡GT810，4G显存哦，只要399。您看那个GTX960
&&&&&&&&& 才2G显存，跑游戏卡死了。
&&好，我就买这个了。
买显卡，看显存大小，似乎已经成了电脑小白必带的标签。
显存是什么？显卡中最慢的外部存储体，焊在PCB板上，与GPU独立，想来几斤就有几斤。
?DDRIII的内存，传输带宽是30~40GB\s，也就是说，每秒能访存30~40GB数据。
?笔记本GDDR5的显存，带宽大概是80GB\s。
?GTX960的带宽是148GB\s。
?淘宝杂牌显卡，带宽大概是40~50GB\s。
?GPU的Cache带宽大概是1.5TB\s，CPU的寄存器带宽大概是15TB\s。
?GPU大概有50%的T周期浪费在访问显存上。作为GPU上最慢最庞大的存储体，显存带宽极为重要。
★★★结论：根据木桶原理，访存速度远远远远远小于并行处理速度，因而I/O瓶颈是GPU的最大问题。
★★★CPU不可以直接访显存，仅__device__、__global__函数有权访问，所以需要内存显存数据相互转移。
3.2.2 内存与显存大挪移
一般来说，一个CUDA程序必然少不了以下三步：
?cudaMalloc：创建新的动态显存堆
?cudaMemcpy：将主机(Host)内存复制到设备(Device)显存
?显存处理完之后，cudaMemcpy：设备(Device)显存复制回主机(Host)内存，释放显存cudaFree
其中第三步最容易遗忘。要知道，CPU最后是无法使用显存中的数据的。
3.2.3 GPU设备选择
CUDA并行有三大境界
?多台机器间并行
?多路显卡间并行
?单卡多线程并行
一般单机中需要操控多路显卡交火，比如AlexNet就是典型的双路交火结构。
不像孤岛危机之类的游戏那样，CUDA需要手动检测、切换设备，这是一切计算的开头准备工作。
下面的InitCUDA函数演示了单卡计算的设备选择：
#include "cuda_runtime.h"
#include "cstdio"
bool InitCUDA()
cudaGetDeviceCount(&cnt);
printf("No GPU Device\n");
return false;
for (i = 0; i & i++)
cudaDeviceP
if (cudaGetDeviceProperties(&device, i)==cudaSuccess)
if (device.major &= 1) break;
if (i == cnt)
printf("GPU's Computing Capability must higher than 1.0\n");
return false;
cudaSetDevice(i);
return true;
InitCUDA()
只有计算能力&=1.0(8800GT以上)的N卡才能做CUDA计算，默认选择符合条件的第一块卡。
3.2.4 HelloWorld
GPU版HelloWorld，主要目的是演示CUDA基本程序框架：
?将HelloWorld复制进显存
?让GPU完成strcpy函数
?将显存中的HelloWorld转回内存，并且打印
/****kernel.cu****/
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
__global__ void cudaStrcpy(char *des, char *src)
/*内核函数*/
while ((*src) != '\0') *des++ = *src++;
*des = '\0';
/****gpu_helloworld.cu****/
#include "device.cu"
#include "kernel.cu"
#include "cstring"
void helloworld(char *str1, char *str2)
InitCUDA();
char *dev_str1=0, *dev_str2=0;
int size = strlen(str1) + 1;
cudaMalloc((void**)&dev_str1, size);
/*cuda系函数必须放在cu文件里*/
cudaMalloc((void**)&dev_str2, size);
cudaMemcpy(dev_str1, str1, size,cudaMemcpyHostToDevice);
cudaStrcpy&&&1,1&&&(dev_str2, dev_str1); /*单线程块、单线程*/
cudaMemcpy(str2, dev_str1, size, cudaMemcpyDeviceToHost);
/****main.cpp****/
#include "cstdio"
#include "cstring"
extern void helloworld(char *str1, char *str2);
int main()
char src[] = "HelloWorld with CUDA";
char *des = new char[strlen(src)+1];
helloworld(src, des);
printf("%s\n", des);
&3.3 GPU版向量加法
向量加法是CUDA 7.0在VS中提供的样例模板，演示了并行算法的经典trick：循环消除。
利用单个线程块中，多个线程并发执行，来消除循环。
时间复杂度估计，不能简单从O(n)迁移到O(1)，因为GPU同时并行量存在限制。
即便是Kepler架构中拥有192SP的SM阵列，理论同时并行量也不过是2048。
/****kernel.cu****/
__global__ void kernel_plus(int *a, int *b, int *c)
int x = threadIdx.x;
c[x] = a[x] + b[x];
/****gpu_vectoradd.cu****/
void vectorAdd(int *a, int *b, int *c,int size)
if (!InitCUDA()) return;
int *dev_a = 0, *dev_b = 0, *dev_c = 0;
cudaMalloc((void**)&dev_a, size*sizeof(int));
cudaMalloc((void**)&dev_b, size*sizeof(int));
cudaMalloc((void**)&dev_c, size*sizeof(int));
cudaMemcpy(dev_a, a, size*sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(dev_b, b, size*sizeof(int), cudaMemcpyHostToDevice);
kernel_plus && &1, size && &(dev_a, dev_b, dev_c);
cudaMemcpy(c, dev_c, size*sizeof(int), cudaMemcpyDeviceToHost);
/****main.cpp****/
extern void vectorAdd(int *a, int *b, int *c, int size);
int main()
int a[5] = { 1, 2, 3, 4, 5 }, b[5] = { 10, 20, 30, 40, 50 }, c[5] = { 0 };
vectorAdd(a, b, c,5);
printf("{1,2,3,4,5} + {10,20,30,40,50} = {%d,%d,%d,%d,%d}\n",
c[0], c[1], c[2], c[3], c[4]);
阅读(...) 评论()当前位置： >>
GPU、CUDA 计算高级优化技术精简手册
GPU编程高级优化技术杂谈前言 数年前，我初入编程领域,一开始根据兴趣GPGPU这个方向, 那是CUDA和OpenCL还未出现，那是底层汇编找色器的时代，而我当时正是通过OpenGL使用GPU汇编指令 ，Cg以及GLSL编写着色器来进行GPU通用计算,直至现在一直从事基于GPU和CPU高性能异构计算的工作。 数年前以网名cyrosly经常混迹于CSDN CUDA论坛和CUDA计算QQ群，讨论各种相关技术或是对一些网友的问题答疑解惑。现在想来那段时期， 有过狂妄，有过激情，更结识了友情。是我的好友郑经维正是在学习CUDA技术的过程中结识的，当 时虽只一面之缘，却成为了这个圈子中最要好的朋友,也正是因为他的劝说才有了我写此书的决定， 或许一本对上感觉的书对于读者的意义远大于在论坛上回答成千上百个问题。本来经纬是想让我尽 可能出版的，但是由于工作的原因，没有多余的精力继续写下去，因此打算把还远未未完成的残稿 贡献出来。本书的目的跳过众多相关书籍频繁重复的内容通过几个有趣的实例直接介绍GPU编程中的高级 优化技术，读者可从本书中一窥诸如cublas,cufft那些高性能库的大概面貌和其中所使用的主要优化技 术。当然,即使是初学者,也可以通过本书达到技术上跳跃式的升级，而作者也信奉一个观念：一看就 懂的书不是好书，因为这代表了读者最终可以从中获取的信息量太少抑或是自己潜意识里早已知晓 但并未显现，并可能引发部分读者的思考：是否物有所值。本书内容绝不雷同，力求精简，节奏很 快，希望读者可以通过分析本书中的代码找到开发高质量GPU程序的感觉。作者本人没有写书的经验 ，，甚至自认不太擅长摆弄文字，所以本书未必是一本好书，但绝对有其独特之处，如果读者可从 中或多或少的学到些其它相关书籍中没有见过的内容，那么也不枉此书了。写作过程历时约一个半 月，由于写作仓促，因此不免有疏漏之处,若有发现,可联系作者更正。作者的联系方式： QQ 微信 : : 目录 第一章 1.0 CUDA设备 1.0.0 核心微架构 1.0.1 寄存器文件结构 设备微架构1.0.2 指令流水线 1.1 GCN设备 1.1.0 核心微架构 1.1.1 寄存器文件结构 1.1.2 指令流水线 1.2 GPU设备上的条件分支 第二章 GPU矩阵乘法的高效实现 2.0 前言 2.1 指令级并行和数据预取 2.2 双缓冲区 2.3 宽数据内存事务 2.4 二级数据预取 2.5 细节调优 第三章 基于GPU的稀疏矩阵直接求解器 3.0 简介 3.1 基于quotient graph的符号分析 3.1.1 顶点重排序 3.1.2 构建消去树 3.1.3 寻找超结点 3.1.4 符号分解 3.2 多波前法 3.3 超节点方法 3.4 多波前+超节点方法的并行分解算法 小结 参考资料 第四章 高性能卷积神经网络的实现 4.0 简介4.1 卷积层的高效计算 4.1.1 基于矩阵乘法的卷积 4.1.2 改进-无需额外存储空间的矩阵乘法卷积 4.1.3 高效的FFT实现 4.1.4 基于FFT的快速卷积 4.2 采样层的高效计算 4.2.1 下采样 4.2.2 上采样 4.3 梯度更新的高效实现 4.3.1 偏置的更新 4.3.2 激活值的更新 第五章 多设备编程建议 第六章 GPU编程优化技术总结 6.1.0 CUDA设备上的优化技术 6.1.1 访存优化 6.1.2 指令优化 6.1.3 内核调用优化 6.2.0 GCN设备上的优化技术 6.2.1 访存优化 6.2.2 指令优化 6.2.3 内核调用优化 6.3 构建性能可移植的程序 小结 参考资料第一章 设备微架构前言第一章我们介绍CUDA设备和GCN设备的微架构做。对设备微架构的了解可以在深度优化时提供 理论依据和方向指导，对微架构方面细节的掌握有时甚至是帮助某些应用达到最优性能必须要 的。当然，对底层架构细节的了解并不是必须的，若读者对这些内容没有兴趣，可以跳过本章 。1.0CUDA设备微架构 kepler架构包含了数组SMX,每个SMX有以下功能单元构成： 1 指令预取单元和微指令缓冲区 2 4个warp调度器,每个warp调度器对应两个指令分派单元 3 192个CUDA Core和64或8个双精度计算单元 4 32个超越函数计算单元(SFU) 5 分支逻辑控制单元 6 32个LD/ST存储器访问单元 7 片上缓存结构，包括共享内存，L1缓存，纹理缓存，常量内存以及只读缓存，不同的设备 大小可能不同kepler设备SMX微架构图各种指令的在不同的功能单元上执行，大致可分为四类：简单计算指令，复杂计算指令，分 支指令和访存指令，下面是各个单元所支持的操作（仅以kepler和maxwell设备为例） CUDA Core ：32位单精度浮点加法，乘法，积和融加运算；32位整数加法；32位数据的比较操作，最小和 最大操作；32位数据的位逻辑操作（and, or, xor）;8位，16位数据和32位数据之间的转换操作。 双精度计算单元：双精度浮点加法，乘法，积和融加运算。 SFU：整数乘法,除法；单精度浮点数除法以及基本的数学函数如正弦，余弦和指数等操作；pop c，clz，brev，bfe和bfi操作。 分支逻辑控制单元：分支，跳转等逻辑操作。 LD/ST单元：全局内存，共享内存，局部内存，常量内存，纹理加载等存储器访问操作。 warp vote和warp shuffle操作时在是专门的组合逻辑单元上完成的。 maxwell设备和kepler类似，但是每个SMM里包含了四组独立的warp计算单元，每个warp单元包含了 一个微指令缓冲区，两个指令分派单元，1个warp调度器，32个CUDA Core，1个双精度计算单元，8个单精度SFU,4个纹理单元，16k个32位寄存器组成的寄存器文件以及24k纹理缓存。所有四个warp单元共享一个指令缓存以及64k~96k的共享/L1缓存。虽然每个SMM中的CUDA Core数量少于SMX，但是每个计算单元具有更高的性能功耗比。相对来说，maxwell具有更优秀的单精度 计算效能，但为了平衡性能功耗比，所有计算能力的maxwell设备对双精度计算的支持都十分有限，而 kepler更适合那些需要双精度计算的专业领域。每三个SMX或每四个SMM组成一个GPC，所有GPC共享512k~ 2M的二级缓存。 缓存小结： 对于计算能力2.x的设备，纹理管线完全独立于L1/L2缓存结构。 对于计算能力3.x的设备，纹理和L1缓存是各自独立的，但都通过L2缓存加载全局内存数据。 对于计算能力5.x的设备，共享内存，纹理缓存都属于L1缓存的一部分，所有对全局数据的访问 都通过L2缓存路径，除非采用直写策略。 maxwell设备微架构图1.0.1寄存器文件结构正因为寄存器索引包含在指令编码中的，且距离计算单元最近，因此其延迟比其它任何类型的内 存都要低，带宽也要远高于其它类型的内存，这也正是为什么充分使用寄存器是某些应用达到峰值 性能所必须且唯一的原因。在计算能力3.5+的kepler设备和maxwell设备上的每个SM中的64k个32bit寄存器文件被划分为4个bank,每个bank中的16k个32bit寄存器共享32个通道(lane)，每个通道32位宽。由于每个线程最多可使用255个寄存器，那么如果使用 单个bank的话，每个源或目标寄存器在指令编码中需要占用8位；如果使用4个bank，那么每个目标寄 存器或是源寄存器只需占用6位，剩下的位数可以留作它用。比如对与四操作数指令（比如FMA）， 如果采用单个bank结构的寄存器文件，那么寄存器索引在指令编码中需要占用32位；而采用4bank结构的寄存器文件的话，寄存器索引在指令编码中只需要占用24位，节省了8位。在32位操作数的 情况下，寄存器文件和4个bank的映射关系为 bank0 : R0, R4, R8 , R12, … bank1 : R1, R5, R9 , R13, … bank2 : R2, R6, R10, R14, … bank3 : R3, R7, R11, R15, … 因此每个源或目标寄存器在指令编码中无需额外的2位表示寄存器的bank，否则将无任何节省。了解 设备的寄存器文件结构对于性能分析以及深度优化具有至关重要的作用，因为不合理的寄存器分配 会造成bank冲突（bankconflict）,类似于共享内存的bank冲突；在没有bank冲突的情况下，指令同时从多个bank存取寄存器数据 。由于源操作数最多为3个，所以寄存器的bank冲突分为1-way,2way两种情况；目标寄存器和元寄存器之间不存在bank冲突的问题，因为取源操作数和写结果操作是 在不同的时钟次序下进行的，。每增加一路冲突，就会增加一个时钟周期的延迟。比如 FMA R1，R7，R3，R1具有1路bank冲突（R7和R3）。如果调整寄存的分配如下 FMA R1，R0，R3，R1则可以消除bank冲突。 使用向量加载可以保证寄存器分配对齐到避免bank冲突的边界上，比如 LD.E.128 R0, [R13] R0，R1，R2，R3四个寄存器分别对应于bank0, bank1, bank2, bank3。注意到R13位于bank1中，但是这并不会和上面代码中的R1产生冲突，因为LD/ST单元会首先取出R 13中的数据作为地址发送到数据总线上，然后（一定的数据传输延迟后）才将数据加载到bank0中的R0 寄存器。使用向量类型时寄存器的分配需要对其到偶数编码的寄存器边界上，因此有时候使用向量 类型会增加寄存器的使用数量，这有点类似于在分配对齐内存时（比如使用AVX256指令时最好将数据 首地址对其到32字节边界上），实际分配的内存空间会比通常的内存分配大一点。 寄存器文件分配需要注意的另一个问题是操作数端口对齐，简单来说就是在为指令分配寄存 器时，相邻的指令中位于同一个bank中的寄存器最好是具有相同的位置，例如FMA R2, R0, R1, R2 FMA R6, R1, R4, R6 {R0,R4}/R1位于同一个bank中，但是在指令中的次序却不同，如果调整为 FMA R2， R0， R1， R2 FMA R6， R4， R1， R6 则可以对齐到操作数端口上，编译器也往往在在指令间插入MOV操作来达到操作数端口对齐的目的 FADD R2, R0, R1 FADD R2, R1, R2 =& FADD R2, R0, R1 MOV R4, R2 FADD R2, R4, R1 在64位操作数的情况下，寄存器和4个bank的映射关系为 bank0 : R[0:1], R[ 8: 9], … bank1 : R[2:3], R[10:11],… bank2 : R[4:5], R[12:13], … bank3 : R[6:7], R[14:15], … 你可能会奇怪，为什么对于32位操作数和64位操作数寄存器和bank之间的映射关系不同。这是因 为bank和寄存器之间并不是隶属关系，bank仅仅相当于数据通道，寄存器文件中的任何数据都可以从 任何一个bank“进出”。可以将bank设想成地铁站的多个出入口，将每个warp中的bank对应的32个通道比作 地铁的各个门，而将人比作数据。那么，人们首先从地铁站的各个出入口进站，然后再从各个地铁 的门进入地铁；反向的情况则是人们从地铁的各个门离开地铁，然后从各个出入口离站。注意，“进 出”的不是寄存器本身，而是寄存器中所包含的数据，寄存器文件本身是一组固定在芯片上的门电路 阵列，是无法移动的；而我们在编程中所说的寄存器指的仅是物理寄存器的标示符；所以为了区分 ，常将两种情况分称为物理寄存器和逻辑寄存器，以下若无特别说明，指的均是逻辑寄存器。 由于每个warp中的32个线程共享四个bank,每个bank有32个32位通道，因此对于64位操作数需要连 续的两个通道或是根据高低32位分两次存取数据，所以每个warp调度器为对应的warp中所有的线程发 射指令需要两个时钟的延迟（或者说每个warp调度器在1个时钟周期只能为16个双精度单元发射指令 ，这样在没有bank冲突的情况下，4个warp调度器在单个时钟周期内可以以全吞吐率为SM内的全部64个 双精度单元进行指令译码）,下面以DFMA操作说明warp中32个线程通过32路通道存取寄存器文件的操作 DFMA R4， R0， R2， R4 第一个时钟周期:lane 0 bank0 bank1 bank212… 31R0, R0, R0, …, R0 R2, R2, R2, …, R2 R4, R4, R4, …, R4第二个时钟周期: lane 0 bank0 bank1 bank2 1 2 … 31R1, R1, R1, …, R1 R3, R3, R3, …, R3 R5, R5, R5, …, R5L1缓存，常量内存，共享内存以及纹理缓存的信息在很多资料中都有详细的叙述，这里不再多说 ，有兴趣的读者可查阅本章后面的参考资料。 1.0.2 指令流水线首先指令预取单元从指令缓存（片上的缓存，物理存储介质和L1缓存相同）中取出多条指令进行 译码，然后将译码后的指令存储在隶属warp调度器的微指令缓冲区中。指令发射单元从微指令缓冲区 中取出两条指令分派给warp调度器选定的warp指令管线中，两条同类指令是顺序进入管线的，指令管 线通过将一条指令的发射，取操作数时钟和另一条指令的取操作数和在功能单元上的执行时钟重叠 来流水线化指令的操作，这样可以有效降低整个指令流水期间所有指令的总时钟数。流水线化在多 个层次都同时存在，如不同功能单元之间的流水化，warp单元和warp单元之间的流水化以及指令和指 令之间的流水化。不同功能单元之间的流水化，比如将计算单元上的操作和LD/ST单元上的操作重叠 进行来隐藏数据访问的延迟。通过多个warp之间的切换让处于停顿延迟中的warp和正在执行的warp的 时钟重叠隐藏warp中的计算或是访存延迟，亦即将warp流水化。单条指令在管线中的过程分为发射， 执行和写回。发射阶段从寄存器文件中取操作数；在执行阶段根据指令中的操作码在相应的单元上 执行特定的操作；写回阶段则将计算结果写入目标寄存器。因此，指令之间的流水化正是通过将前 一条指令的执行和写回与下一条指令的发射和执行的时钟重叠来减少整个管线中所有指令的操作延 迟。在指令译码阶段完成数据的存取后指令的计算结果并非立即可用，而是需要数个时钟的延迟。 一般来说指令管线的深度至少不应小于指令的计算延迟，才足以保证指令流过管线后已经完成计算 。理想情况下，每当一条指令从指令管线末端流出（完成计算）时，管线的前端就会流入一条新的 指令，这样指令管线中的多条指令通过时钟重叠的方式有效的隐藏计算延迟。不同的指令可能具有 不同的指令发射周期，亦即当发射一条指令后，需要几个时钟的停顿后才会发射下一条同类指令。 我们假设一个指令管线的深度和指令的执行延迟均为4个时钟周期，每周期均可发射1条指令，那么 指令流水线情况如下 1 clock I0 2 clock I1 I0 3 clock I2 I1 I0 4 clock I3 I2 I1 I05 clock 6 clock 7 clock 8 clockI3 I2 I1 I0 I3 I2 I1 I0 I3 I2 I1 I0 I3 I2 I1 I0 流水线化,红色部分是从管线末端流出的以完成操作如果没有采用流水线化，那么 1~ 5 clock I0 6~10 clock 11~15 clock 16~20 clock I1 I2 I3可以看出，采用流水线化的操作完成四条指令的译码和执行只需要8个时钟周期，而不采用流水线化 的操作则需要20个时钟周期。 1.1 GCN设备微架构GCN设备由多组Compute Unit(简称CU),每组CU中包含组4个16路向量SIMD单元，一个标量单元和一个分支通信单元。每个SIMD单元 配备64K的向量寄存器（寄存器粒度是32位），每个线程最多可以使用255个寄存器；向量单元主要执 行简单浮点和整数计算，如32位和64位的浮点加法，浮点乘法,FMA以及32位整数加减法以及24位整数乘 法和乘加操作。标量单元配备8K的标量寄存器以及一个整数计算单元，寄存器的粒度同样是32位大小 。标量单元负责流控制等操作，而整数单元则负责全精度的整数运算。每个CU中的wavefront调度器每 个时钟周期可以发射5条指令，包括4条向量指令（每个SIMD单元一条）和一条标量指令。多个SIMD单 元可以运行同一个workgroup中的不同线程，但是同一个wavefront中的64个线程智能在同一个SIMD单元上分4个时钟周期执行， 即使每个workgroup中只有一个wavefront,线程调度器也不会将线程平均分配到CU中的各个SIMD单元上：GCN设备上的最 小并发调度粒度是wavefront(目前是64)。关于GCN架构更详细的说明可以参考AMD的GCN架构白皮书。 1.1.1 寄存器文件结构 GCN设备上的每个Compute Unit具有256k个32bit寄存器，寄存器文件被划分成4个bank,各自分属Compute Unit内的4组向量计算单元使用。每个bank有64（不同的设备可能不同）个通道，通道宽度是4字节，一 个向量计算单元中每个ALU最多可以使用4个通道（4操作数指令）。在GCN设备上没有寄存器bank冲突 的问题，因为4个bank各自独立，每个向量计算单元只能使用与其对应的一个bank。当使用8字节或16字 节的数据时，使用连续的2个或4个寄存器，但是寄存器的分配要满足对其要求。对于8字节数据，寄 存器编号需要对齐到2的倍数；对于16字节数据，寄存器编号需要对齐到4的倍数；因此当使用向量类 型时有时会增加寄存器的使用数量，这一点类似于对齐的内存分配会比普通的分配略大（比如使用A VX256指令时，通常要求存取数据的首地址按照32字节对齐）。 1.1.2 指令流水线 GCN设备上的指令流水线原理类似于CUDA设备，因此这里不再重复。1.2 GPU设备上的条件分支 这里并不假定任何型号的设备，因此这一节所讲的同时适用于CUDA设备，GCN设备和其它宽向 量设备。每个warp或wavefront拥有一组条件掩码寄存器堆栈（堆栈中包含了多个掩码寄存器，每个掩 码寄存器32位（CUDA）或64位（GCN），每位对应一个线程的掩码），GCN设备上的掩码寄存器堆栈深 度为6。假设我们有一个深度为为3的嵌套分支，那么warp或wavefront中线程的执行过程为 1 首先将计算得到的最外层的条件掩码写入掩码寄存器堆栈顶部的第1个寄存器(掩码堆栈指针当 前指向的位置)。 2 掩码寄存器堆栈指针减1，接着将计算得到的第2层分支的条件掩码写入掩码寄存器堆栈的下 一个寄存器(掩码堆栈指针当前指向的)。 3 掩码寄存器堆栈指针减1，接着将计算得到的第3层分支的条件掩码写入掩码寄存器堆栈的下 一个寄存器。 4 从掩码寄存器堆栈指针指向的位置取出掩码，如果warp或wavefront中的线程号对应的位值为1， 则执行第3层分支内的计算，否则忽略。 5 掩码寄存器堆栈指针加1，取出掩码寄存器中的掩码，如果warp或wavefront中的线程号对应的位 值为1，则执行第2层分支内的计算，否则忽略。6 掩码寄存器堆栈指针加1，取出掩码寄存器中的掩码，如果warp或wavefront中的线程号对应的位 值为1，则执行第1层分支内的计算，否则忽略。 逻辑分支图小结 本章我们分别介绍了CUDA设备和GCN设备的微架构，重点分析了寄存器文件结构和指令管线以 及GPU设备上实现条件分支的原理。深入了解这些内容不仅可以让你知道指令在硬件中具体的行为以 及为什么会这样，同时还会在优化过程中的某个角落为你说明性能无法进一步提升的障碍所在。通 常理解和使用所有这些内容需要一个漫长的过程，建议经常到设备商家官方网站上下载最新的技术 资料，同时一些编译原理和计算机结构方面的书籍也会帮到你，因为很多技术原理无论在CPU还是GPU 上都是通用的。 参考资料 1《CUDA Programming Guide》 2《CUDA Binary Utilities》 3《NVIDIA Kepler GK110 Whitepaper》 4《NVIDIA Geforce GTX980 Whitepaper》 5《AMD GCN Architecture Whitepaper》 6《AMD Southern Island Series Instruction Set Architecture》 7《Graphics Core Next Architecture，Generation 3》 8《Compute Systems A programmer’s Perspective》，Randal E.Bryant, David R.O’Hallaron 9 Performance Upper Bound Analysis and Optimization of SGEMM on Fermi and Kepler GPUs，Junjie Lai, Andr_e Seznec第二章 GPU矩阵乘法的高效实现前言 本章通过介绍开发GPU上的高效矩阵乘法的各种优化方法来展示一些GPU高级优化技术，其中所 用到的大多数优化技术不仅在支持CUDA的设备上有效，对于GCN设备,MIC设备和CPU上同样有效。如非 特别说明，在所有版本的实现中，数据都是列序存储,。2.0 指令级并行和数据预取 对于类似矩阵乘法这样的高计算密度类型的应用，努力提高每个线程中的ILP(instruction level parallism,指令级并行)比单纯使TLP(thread level parallism,线程级并行)更加有效,而且是达到接近峰值性能唯一的途径。ILP的真髓在于通过充分发掘单 一线程内的指令流水线中的并行和并发性（如果看过第一章就会明白，单一线程内的并行存在于 被双发的指令在不同功能单元上的操作，比如在LD/ST上的访存操作和计算单元上的FMA操作可以同 时进行；而并发则存在于相同功能单元上的不同阶段中分时时钟序列中）以及利用寄存器的最低 延迟和超高带宽的优势（在kepler和GCN设备上，寄存器的带宽大约是共享内存的6倍）。循环展开通 常是提高ILP的有效方法，原因在于通过循环展开编译器可以在更长的无数据依赖的指令序列中对 指令进行排序从而最大化指令吞吐率并最小化流水线中的停顿，同时减少循环索引和部分存储器 地址的计算。 dgemm, loop unrolling&prefetch __device__ __forceinline__ void d_rank8x8( double* C, const double* A, const double* B ) { double a[8], a[0]=A[0*16]; a[1]=A[1*16]; a[2]=A[2*16]; a[3]=A[3*16]; a[4]=A[4*16]; a[5]=A[5*16]; a[6]=A[6*16]; a[7]=A[7*16]; #pragma unroll for( int i=0; i&8; ++i ){ b=B[i*16]; C[i*8+0]+=a[0]*b; C[i*8+1]+=a[1]*b; C[i*8+2]+=a[2]*b; C[i*8+3]+=a[3]*b; C[i*8+4]+=a[4]*b; C[i*8+5]+=a[5]*b; C[i*8+6]+=a[6]*b; C[i*8+7]+=a[7]*b; } } __global__ void cuk_dgemm_unroll( double* d_C, const double* d_A, const double* __restrict__ d_B, int n, int lda, int ldb, int ldc ) { __shared__ double smem[2048]; double p0, p1, p2, p3, q0, q1, q2, q3, c[64]={0.f}; int k, lane, lane=threadIdx.x&15; slot=(threadIdx.y&&1)+(threadIdx.x&&4); d_C+=(((blockIdx.y&&7)+slot)*ldc+(blockIdx.x&&7)+lane); d_A+=(threadIdx.y*lda+(blockIdx.x&&7)+threadIdx.x); d_B+=(((blockIdx.y&&7)+((threadIdx.x&1)&&4)+(threadIdx.x&&1))*ldb+threadIdx.y);double* St=&smem[(threadIdx.y&&7)+threadIdx.x]; double* At=&smem[lane]; double* Bt=&smem[1024+slot]; if(threadIdx.y&n) { p0=d_A[0*32]; p1=d_A[1*32]; p2=d_A[2*32]; p3=d_A[3*32]; q0=d_B[0*32*ldb]; q1=d_B[1*32*ldb]; q2=d_B[2*32*ldb]; q3=d_B[3*32*ldb]; } for( k=n-8; k&=0; k-=8 ) { *(St+0*32)=p0; *(St+1*32)=p1; *(St+2*32)=p2; *(St+3*32)=p3; *(St+0*32)=q0; *(St+1*32)=q1; *(St+2*32)=q2; *(St+3*32)=q3; __syncthreads(); if(threadIdx.y&k){ d_A+=(lda&&3); d_B+=8; p0=d_A[0*32]; p1=d_A[1*32]; p2=d_A[2*32]; p3=d_A[3*32]; q0=d_B[0*32*ldb]; q1=d_B[1*32*ldb]; q2=d_B[2*32*ldb]; q3=d_B[3*32*ldb]; } #pragma unroll for( int i=0; i&8; ++i ){ d_rank8x8( c, &At[i*128], &Bt[i*128] ); } __syncthreads(); } if(k!=-8) { *(St+0*32)=p0; *(St+1*32)=p1; *(St+2*32)=p2; *(St+3*32)=p3; *(St+0*32)=q0; *(St+1*32)=q1; *(St+2*32)=q2; *(St+3*32)=q3; __syncthreads(); do{d_rank8x8(c, At, Bt); At+=128; Bt+=128; }while((++k)&=0); } #pragma unroll for( int i=0; i&64; i+=8 ){ d_C[0*16]=c[i+0]; d_C[1*16]=c[i+1]; d_C[2*16]=c[i+2]; d_C[3*16]=c[i+3]; d_C[4*16]=c[i+4]; d_C[5*16]=c[i+5]; d_C[6*16]=c[i+6]; d_C[7*16]=c[i+7]; d_C+=(ldc&&4); } } 第一个版本中我们主要的优化方法是通过循环展开获得更多的ILP以及通过数据预取和计算的 重叠隐藏访存延迟。（注：循环展开并非在任何情况下都能提升指令级并行，待展开的计算中还需 要尽可能的无数据依赖关系，待展开的循环体中的计算也应尽可能的简单；而且过度的循环展开还 可能会造成寄存器溢出或是设备资源占用率过低从而导致性能下降）。 注意内核参数“const double* __restrict__ d_B”,对于计算能力3.5+的设备，编译器会自动使用LDG指令通过纹理管线加载指令加载数据，从而可以 利用纹理缓存较少非合并访问的影响。对于计算能力2.x或3.x的设备,显式的使用纹理通常比普通的缓 存效果更好，虽然很轻微(受限于L2缓存有限的带宽和纹理操作的长延迟)。这是因为纹理管线独立于 普通的存储路径,可以提升可用的缓存带宽(纹理缓存的带宽+L1缓存的带宽)，从而更加平衡负载,并减 小不同数据之间的L1缓存污染(通过纹理缓存加载的数据不经过L1缓存)。对于计算能力5.x的设备，由 于统一了L1缓存和纹理缓存，因此使用线性纹理和普通L1缓存已经没有本质的区别了。 此外所选择的的网格布局也会直接影响性能和代码的可读性。根据不同的架构选择每线程计 算的元素数量，太少的话不能充分利用寄存器文件高带宽低延迟的优势；太多的话会造成寄存器溢 出。对于计算能力2.x和3.0的设备，每个线程计算16个元素(1x16或者4x4),线程块的大小选择为256(计算能 力2.x，block布局可设置为64x4或32x8,每个block计算64x64的子矩阵)或1024(计算能力3.0，block布局可设置为128 x8或32x32,每个block计算128x128的子矩阵)。对于计算能力3.5+,5.x的设备，比较好的选择是每线程计算8x8 个元素(对于比较小的规模可以每个线程计算4x8个元素，每个block计算128x64大小的子矩阵)，线程块的 大小设置为256(block布局可设置为16x16,32x8,128x2或者256；每个block计算128x128的子矩阵,一个SM最多可支 持一个block并发(对于计算能力3.7的设备,每个SMX上的活跃block队列最多可部署2个)。 我们采用32x8的block布局，每个线程计算64个元素，每个block计算128x128的子矩阵。tx ty0 ty1 ty2 ty3 ty4 ty5 ty6 tx + 32 tx + 64 tx+96ty0 ty1 ty2 ty3 ty4 ty5 ty6 ty7s0 lane 16 + lane 32 + lane 48 + lane 64 + lane 80 + lane 96 + lane 112 + lanes1s2s3s4s5s6s7s8s9s10 s11s12 s13 s14s152.1 双缓冲区技术 我们知道GPU网格调度器通过block内的warp(在AMD的GPU中称之为wavefront)切换隐藏内存访问和计 算延迟；而通过block之间的切换隐藏块内同步延迟，但这里由于寄存器和共享内存资源的限制，一 个SMX上最多只能部署一个活跃block（对于计算能力3.7的设备，由于具有双倍的共享内存和寄存器文 件，所以每个SMX可以支持两个block并发），因而无法再通过增加每个SMX上可以并发的block数量来进 一步提升效率。类似单个线程内的并行和并发，这时我们就需要考虑如何在单个block内寻找更多的 并行和并发机制。容易看出主循环体内的并行和并发受限于块内同步，因此如果消除下一个循环步 预取的数据和当前计算中所用数据的依赖关系，则第二个同步原语是可以消除的，为此我们通过双 倍的共享内存作为双缓冲来消除这种影响： dgemm,loop unrolling+double buffering __global__ void cuk_dgemm_unroll_db( double* d_C, const double* d_A, const double* __restrict__ d_B, int n, int lda, int ldb, int ldc ) { __shared__ double smem[4096]; … if(threadIdx.y&n) { smem[0*32+sidx]=d_A[0*32]; smem[1*32+sidx]=d_A[1*32]; smem[2*32+sidx]=d_A[2*32]; smem[3*32+sidx]=d_A[3*32]; smem[+sidx]=d_B[0*32*ldb]; smem[+sidx]=d_B[1*32*ldb]; smem[+sidx]=d_B[2*32*ldb]; smem[+sidx]=d_B[3*32*ldb]; } for( k=n-8; k&=0; k-=8 ) {__syncthreads(); if(threadIdx.y&k){ d_A+=(lda&&3); d_B+=8; sidx^=0x800; smem[0*32+sidx]=d_A[0*32]; smem[1*32+sidx]=d_A[1*32]; smem[2*32+sidx]=d_A[2*32]; smem[3*32+sidx]=d_A[3*32]; smem[+sidx]=d_B[0*32*ldb]; smem[+sidx]=d_B[1*32*ldb]; smem[+sidx]=d_B[2*32*ldb]; smem[+sidx]=d_B[3*32*ldb]; } #pragma unroll for( int i=0; i&8; ++i ){ d_rank8x8( c, &smem[i*128+lane], &smem[i*128+slot] ); } lane^=0x800; slot^=0x800; } …… } 使用双缓冲技术将当前计算的读取缓冲区和负责装载下一步数据的缓冲区相分离，从而加载 下一步数据到共享内存的操作可以和当前的计算重叠执行，减少了同步开销从而减少流水线的中断 。但是使用双缓冲时需要注意当前的设备上具体的效果，此外不同的计算规模效果也会不同。对于 计算能力2.x的设备，由于双精度指令和存储指令不能双发(因此LD/ST指令和双精度浮点计算指令是顺 序发射的，但在指令完成操作的这段延迟期间，其执行仍是重叠的)，在其上使用双缓冲区的效果相 比kepler和maxwell设备会小些。 2.2 宽内存事务 我们还可以做得更好,对内层循环(d_rank8x8)进行分析可知，每个循环步共有512条FMA指令和128条 共享内存加载指令，浮点计算指令所占比例(这里我们不考虑数据预取和缓冲区地址计算的相关指令, 这样并不影响分析的有效性)是0.8。如果能够进一步提升这个比例，就有可能进一步提升性能。因为F MA操作的数量是固定的，因此可以通过使用128bit存储操作来缩减每个线程的内存事务的数量。首先 需要对d_rank8x8进行改写： #define mRank8x8(i0,i1,i2,i3){ \ c[0*8+0]+=a[i0].x*b[i0].x;\ c[0*8+1]+=a[i0].y*b[i0].x;\ c[0*8+2]+=a[i1].x*b[i0].x;\ c[0*8+3]+=a[i1].y*b[i0].x;\ c[0*8+4]+=a[i2].x*b[i0].x;\ c[0*8+5]+=a[i2].y*b[i0].x;\ c[0*8+6]+=a[i3].x*b[i0].x;\ c[0*8+7]+=a[i3].y*b[i0].x;\ c[1*8+0]+=a[i0].x*b[i0].y;\ c[1*8+1]+=a[i0].y*b[i0].y;\ c[1*8+2]+=a[i1].x*b[i0].y;\ c[1*8+3]+=a[i1].y*b[i0].y;\ c[1*8+4]+=a[i2].x*b[i0].y;\ c[1*8+5]+=a[i2].y*b[i0].y;\ c[1*8+6]+=a[i3].x*b[i0].y;\c[1*8+7]+=a[i3].y*b[i0].y;\ c[2*8+0]+=a[i0].x*b[i1].x;\ c[2*8+1]+=a[i0].y*b[i1].x;\ c[2*8+2]+=a[i1].x*b[i1].x;\ c[2*8+3]+=a[i1].y*b[i1].x;\ c[2*8+4]+=a[i2].x*b[i1].x;\ c[2*8+5]+=a[i2].y*b[i1].x;\ c[2*8+6]+=a[i3].x*b[i1].x;\ c[2*8+7]+=a[i3].y*b[i1].x;\ c[3*8+0]+=a[i0].x*b[i1].y;\ c[3*8+1]+=a[i0].y*b[i1].y;\ c[3*8+2]+=a[i1].x*b[i1].y;\ c[3*8+3]+=a[i1].y*b[i1].y;\ c[3*8+4]+=a[i2].x*b[i1].y;\ c[3*8+5]+=a[i2].y*b[i1].y;\ c[3*8+6]+=a[i3].x*b[i1].y;\ c[3*8+7]+=a[i3].y*b[i1].y;\ c[4*8+0]+=a[i0].x*b[i2].x;\ c[4*8+1]+=a[i0].y*b[i2].x;\ c[4*8+2]+=a[i1].x*b[i2].x;\ c[4*8+3]+=a[i1].y*b[i2].x;\ c[4*8+4]+=a[i2].x*b[i2].x;\ c[4*8+5]+=a[i2].y*b[i2].x;\ c[4*8+6]+=a[i3].x*b[i2].x;\ c[4*8+7]+=a[i3].y*b[i2].x;\ c[5*8+0]+=a[i0].x*b[i2].y;\ c[5*8+1]+=a[i0].y*b[i2].y;\ c[5*8+2]+=a[i1].x*b[i2].y;\ c[5*8+3]+=a[i1].y*b[i2].y;\ c[5*8+4]+=a[i2].x*b[i2].y;\ c[5*8+5]+=a[i2].y*b[i2].y;\ c[5*8+6]+=a[i3].x*b[i2].y;\ c[5*8+7]+=a[i3].y*b[i2].y;\ c[6*8+0]+=a[i0].x*b[i3].x;\ c[6*8+1]+=a[i0].y*b[i3].x;\ c[6*8+2]+=a[i1].x*b[i3].x;\ c[6*8+3]+=a[i1].y*b[i3].x;\ c[6*8+4]+=a[i2].x*b[i3].x;\ c[6*8+5]+=a[i2].y*b[i3].x;\ c[6*8+6]+=a[i3].x*b[i3].x;\ c[6*8+7]+=a[i3].y*b[i3].x;\ c[7*8+0]+=a[i0].x*b[i3].y;\ c[7*8+1]+=a[i0].y*b[i3].y;\ c[7*8+2]+=a[i1].x*b[i3].y;\ c[7*8+3]+=a[i1].y*b[i3].y;\ c[7*8+4]+=a[i2].x*b[i3].y;\ c[7*8+5]+=a[i2].y*b[i3].y;\ c[7*8+6]+=a[i3].x*b[i3].y;\ c[7*8+7]+=a[i3].y*b[i3].y;\ } #define mFetchSmem(i0,i1,i2,i3,k){ a[i0]=smem[k*64+ 0+lane]; a[i1]=smem[k*64+16+lane]; \ \ \a[i2]=smem[k*64+32+lane]; a[i3]=smem[k*64+48+lane]; b[i0]=smem[512+k*64+ 0+slot]; b[i1]=smem[512+k*64+16+slot]; b[i2]=smem[512+k*64+32+slot]; b[i3]=smem[512+k*64+48+slot]; }\ \ \ \ \ \通过使用128bit存储器操作，我们将共享内存加载指令数目减少了一半，现在只有64条共享内 存加载指令，浮点计算指令所占比例约是0.89。虽然128bit的共享内存操作以及其引发的bank conflicts会带来更高的延迟，但是通过合理的安排计算足以将这些延迟抵消掉。 dgemm, loop unrolling+double buffering+128bit LD/ST __global__ void cuk_dgemm_unroll_db_128b( double* d_C, const double* d_A, const double* __restrict__ d_B, int n, int lda, int ldb, int ldc ) { __shared__ double2 smem[2048]; double2 a[4], b[4]; …… if(threadIdx.y&n){ ((double*)smem)[0*32+sidx]=d_A[0*32]; ((double*)smem)[1*32+sidx]=d_A[1*32]; ((double*)smem)[2*32+sidx]=d_A[2*32]; ((double*)smem)[3*32+sidx]=d_A[3*32]; ((double*)smem)[+sidx]=d_B[0*32*ldb]; ((double*)smem)[+sidx]=d_B[1*32*ldb]; ((double*)smem)[+sidx]=d_B[2*32*ldb]; ((double*)smem)[+sidx]=d_B[3*32*ldb]; } for( k=n-8; k&=0; k-=8 ) { __syncthreads(); if(threadIdx.y&k){ d_A+=(lda&&3); d_B+=8; sidx^=0x800; ((double*)smem)[0*32+sidx]=d_A[0*32]; ((double*)smem)[1*32+sidx]=d_A[1*32]; ((double*)smem)[2*32+sidx]=d_A[2*32]; ((double*)smem)[3*32+sidx]=d_A[3*32]; ((double*)smem)[+sidx]=d_B[0*32*ldb]; ((double*)smem)[+sidx]=d_B[1*32*ldb]; ((double*)smem)[+sidx]=d_B[2*32*ldb]; ((double*)smem)[+sidx]=d_B[3*32*ldb]; } #pragma unroll for( int i=0; i&8; ++i ){ mFetchSmem(0,1,2,3,i) mRank8x8(0,1,2,3) } lane^=0x400; slot^=0x400; } ……} 合理使用128bit存储操作可以减少单个线程内的访存事务的数量(也会因此提升计算指令的整体 比例)从而提升性能。对于双字宽数据,提升效果相对较小;因为相对于单字宽数据,双字宽内存事务的 数量减少相对较低,因此计算指令所占比例的提升也较小(不过对于矩阵乘法，这仍然是进一步提升效 率的有效且必须的技术)，在单精度的情况下会得到更好的效果。 2.3 二级数据预取 除了对全局内存的数据进行预取外，我们可以进一步加入共享内存的预取流水,从而通过将计 算和共享内存操作重叠来进一步隐藏访问共享内存的延迟，同时减少计算中等待数据的流水线停顿 。另外，对于双精度矩阵乘法在某些设备上使用双缓冲反而会降低性能,尤其是规模比较小的矩阵不 建议使用；而加入对共享内存数据的预取在多数设备和规模上均能性能提升(作者在GTX680，GTX780，G TX970上测试中都获得了提升)。 dgemm, loop unrolling + double buffering + 128bit LD/ST + smem_prefetch __global__ void cuk_dgemm_unroll_db_128b_prefsmem( double* d_C, const double* d_A, const double* __restrict__ d_B, int n, int lda, int ldb, int ldc ) { …… ((double*)smem)[0*32+sidx]=d_A[0*32]; ((double*)smem)[1*32+sidx]=d_A[1*32]; ((double*)smem)[2*32+sidx]=d_A[2*32]; ((double*)smem)[3*32+sidx]=d_A[3*32]; ((double*)smem)[+sidx]=d_B[0*32*ldb]; ((double*)smem)[+sidx]=d_B[1*32*ldb]; ((double*)smem)[+sidx]=d_B[2*32*ldb]; ((double*)smem)[+sidx]=d_B[3*32*ldb]; __syncthreads(); mFetchSmem(0,2,4,6,0) for( k=n-8; k&=0; k-=8 ) { if(threadIdx.y&k) { d_A+=(lda&&3); d_B+=8; sidx^=0x800; ((double*)smem)[0*32+sidx]=d_A[0*32]; ((double*)smem)[1*32+sidx]=d_A[1*32]; ((double*)smem)[2*32+sidx]=d_A[2*32]; ((double*)smem)[3*32+sidx]=d_A[3*32]; ((double*)smem)[+sidx]=d_B[0*32*ldb]; ((double*)smem)[+sidx]=d_B[1*32*ldb]; ((double*)smem)[+sidx]=d_B[2*32*ldb]; ((double*)smem)[+sidx]=d_B[3*32*ldb]; } mFetchSmem(1,3,5,7,1) mRank8x8(0,2,4,6) mFetchSmem(0,2,4,6,2) mRank8x8(1,3,5,7) mFetchSmem(1,3,5,7,3)mRank8x8(0,2,4,6) mFetchSmem(0,2,4,6,4) mRank8x8(1,3,5,7) mFetchSmem(1,3,5,7,5) mRank8x8(0,2,4,6) mFetchSmem(0,2,4,6,6) mRank8x8(1,3,5,7) mFetchSmem(1,3,5,7,7) mRank8x8(0,2,4,6) __syncthreads(); lane^=0x400, slot^=0x400; mFetchSmem(0,2,4,6,0) mRank8x8(1,3,5,7) } …… } 2.4 细节调优 以上代码都假设A的行数，B的列数以及C的行列数是128的倍数。所以如果矩阵大小不匹配block 的计算尺寸，那么更好的实现是添加一个边界处理的内核专门处理边界。如果边界大小和block的计 算尺寸相差不大(因此不会浪费太多显存空间),则可以通过填充零将矩阵补齐到匹配block的计算尺寸,这 样就不再需要单独的边界处理内核，效率也会更高些。还有一点需要说明：将数据存取操作和浮点 计算操作混合排列可以获得更高的效率（更高的IPC）,这里为了更清晰的表达所以没有这样做。上述 代码的性能仍有较大的提升空间，如果想获得接近峰值的效率，需要手工使用PTX或SASS对细节进行 调优，否则仅仅使用高层的C/C++语言进行内核开发是无法获得接近于峰值性能的。如果开发人员 想通过使用PTX显式的控制诸如寄存器分配，指令排序等问题，可能得不到想要的结果。原因是PTX并 非真正的面相本地机器指令的汇编语言，而是一种为了兼容不同计算能力的设备而设计的一种类似 寄存器传输语言（RTL）的中间层伪汇编语言；从源代码到PTX的编译过程只会做一些初级的优化，多 数优化，比如寄存器分配，指令重排，基本块融合等都是在从PTX代码到SASS原生汇编代码的过程中 进行的(PTX和SASS的关系类似AMD GPU编程中的IL和GCN原生ISA的关系，这里所说的多数情况在IL和GCN ISA之间也同样存在)。因此若要显式的控制指令的执行顺序和寄存器分配等就需要直接使用SASS,不过N VIDIA现在并未给出直接使用汇编编写GPU内核程序的开发工具，也未开放任何关于指令集编码信息的 文档，因此开发者需要自己想办法绕过这个限制。至少对于矩阵乘法，ptxas做的并不是太好（使用C/ C++只能得到最高约75%峰值的效率，使用PTX可以超过80%的峰值，但是离最好的效果还有差距）， 主要有两个原因： 1 指令的排序效果不理想，也就是LD/ST指令，FMA，IADD，LOP， ISETP等指令之间的排列位置不够理想，最好的指令之间的排列距离要综合考虑指令管线的深 度，dual issue,以及各种操作的发射延迟和计算延迟；对指令放置在不同位置进行测试选优也往往是必 不可少的。对于这些情况人工通常可以比编译器做的更好。 2 寄存器bank冲突(bankconflict)和功能单元的操作数端口对齐问题对性能的影响也是非常大的。寄存器bank冲突分为2way和3-way两种情况，对于2-way,每条FMA指令发射到流水线会增加一个时钟的延迟，3way则会增加2个时钟周期，这可能会严重影响效率（每个循环步增加的时钟周期=存在端口冲 突的指令数 x1或是2）,ptxas为了避免寄存器bank冲突以及强制将寄存器对齐到指令的操作数端口而增加了很 多冗余的MOV指令（循环体的追尾效应，使用cuobjdump查看SASS代码会看到在循环体的末尾处多 出很多MOV操作）,虽然MOV操作是直接在寄存器文件内部移动数据（意味着无需像其它操作那 样先从数据端口读取寄存器中的数据，然后再将数据通过数据端口写入寄存器），但仍然不 是免费的,哪怕其延迟很低。即使这样仍然不能保证完全消除寄存的端口冲突。循环体的追尾 效应虽然可以通过在循环体前端移除数据预取代码中的条件代码以及通过循环剥离移除循环 体内的条件代码来消除；但由于第一因素仍然存在，所以性能的提升仍然受限。 最后需要说明的是由于每个线程的输出是8字长，所以每16个线程输出一列不会有非合并写出的问题( 对于同一个warp对应的8字节数据，会分两个批次进行传输，每次传输连续的128字节序列，对应16个 线程)，因此是否通过共享内存对数据进行重新打包对性能影响不大，因此我们并没有使用，但是当 使用单精度实现时需要这样做。以下是使用共享内存对数据进行重打包的示例代码（d_C的地址计算 也需要做相应的改变）： double2* Cs=&smem[(slot&&4)+lane]; volatile double* Ct=&((double*)smem)[(threadIdx.y&&8)+threadIdx.x]; #pragma unroll for( int i=0; i&64; i+=8 ) { Cs[0*16]=make_double2(c[i+0],c[i+1]); Cs[1*16]=make_double2(c[i+2],c[i+3]); Cs[2*16]=make_double2(c[i+4],c[i+5]); Cs[3*16]=make_double2(c[i+6],c[i+7]); d_C[0*32]=Ct[0*32]; d_C[1*32]=Ct[1*32]; d_C[2*32]=Ct[2*32]; d_C[3*32]=Ct[3*32]; d_C[ldc+0*32]=Ct[4*32]; d_C[ldc+1*32]=Ct[5*32]; d_C[ldc+2*32]=Ct[6*32]; d_C[ldc+3*32]=Ct[7*32]; d_C+=(ldc&&4); }下面给出各个版本的测试结果(每个测例循环100次取平均，单位ms) version size cublas unroll unroll+128b unroll+128b+prefsmem unroll+db+128b 0.017000 GTX970 48x96 0.................042550unroll+db+128b+prefsmem本章小结 通过合理组合使用循环展开，128bit存储操作以及对全局内存和共享内存的两级数据预取可以 有效提升矩阵乘法的效率。在这些优化技术中，循环展开，数据预取通常是有效的，而128bit内存操 作和双缓冲技术则需视情况而定。此外，对一些细节的处理也必不可少，比如将固定的索引计算移 除循环并加到指针变量上以减少地址计算和指令数量。当矩阵规模比较小时(比如512x512)，可以每个 线程计算32个元素，每个block计算128x64的输出，从而充分利用硬件资源平衡计算负载。当然，这些优 化技术并不仅限于矩阵乘法。小结 一开始我们就从一个性能很高的版本开始更进一步的优化，之所以没有从一个更加简单的版 本(比如CUDA programming guide上的例子)开始逐步走完全程，并不是因为作者太懒(好吧，我承认是有这个原因)。其一是因为各 种相关书籍都有介绍，除了占用篇幅，作者实在找不到把那些内容放在本书的理由。更主要的原因 是作者更建议一开始就从尽可能高效的方法去实现，这个建议对任何项目都适用。假设你正在做一 个比较大的项目，可能会想着先让项目运行起来，至于效率以后可以逐步改进。这种想法对有些项 目可能适用。但是对于更多的项目这样做可能会大幅度拖延开发进度，因为当你准备优化时，可能 会发现为了改进性能，不得不改变数据结构 (因为有时候大幅度的性能提升必须依赖特定的数据结构和算法)和大范围的改写代码逻辑。有些情况 下一个简单的数据转换接口即可在时间和效率上得到折中（通常意味着不是最好)；但很多时候对数 据结构或是算法的改变却会牵一发而动全身甚至从架构到代码都需要重新设计开发。所以强烈建议 从一开始就用尽可能优化的方法作为起点，用尽可能优质的实现构建你的项目，这样前期可能会需 要多些时间，但是之后会发现这样可以大大提高整体的开发进度，同时项目的质量也会得到保障。 总之，对一个低质量的程序通过渐进的方法不断优化的工作效率远低于直接使用优化的方法；正所 谓快刀斩乱麻，有时候屏蔽已有的东西，可以避免被潜意识中经验性或现有的代码规则所限制。参考资料 1《cuda programming guide》 2《kepler tuning guide》 3《maxwell tuning guide》 4《PTX ISA version 4.2》 5《nvidia-kepler-gk110-architecture-whitepaper》 6 Performance Upper Bound Analysis and Optimization of SGEMM on Fermi and Kepler GPUs，by Junjie Lai, Andr_e Seznec。 7 Anatomy of High-Performance Matrix Multiplication，KAZUSHIGE GOTO,The University of Texas at Austin and ROBERT A. VAN DE GEIJN The University of Texas at Austin。 8 On Reducing TLB misses in Matrix Multiplication,kazushige goto, Robert van de geijn, november 1, 2002。第三章 基于GPU的稀疏直接求解器前言 本章可能是所有章节中最难得了，如果节奏过快，对于即使有一定相关开发经验的人来说要 彻底理解其中的内容也会相当困难，所以作者会给出相比其它章节更多的说明。第一章中类似的技 术将会在本章得到应用，通过阅读本章，读者不仅可以进一步巩固第一章中学到的GPU计算的优化技 术，还能了解到一些有趣的算法。在这个过程中读者还会更深的体会到从算法的伪代码到实际的可 执行代码之间的难度。3.0 稀疏矩阵直接解法概述 很多科学计算问题中都涉及到求解线性方程组，而这些方程组不仅规模很大且通常是稀疏的 。解决这类问题通常使用迭代解法或直接解法。迭代解法一般内存空间需求小，程序实现比较简单 ；但是迭代解法的收敛速度依赖于矩阵的性态，如果矩阵过于病态，即使采用复杂的预处理技术(预 处理越复杂，也意味着求解效率越低，内存空间需求越大，从而迭代解法优势也会逐渐丧失)也很难 收敛甚至不收敛。直接解法的缺点是内存空间的需求比迭代解法大很多，但是其具有解的精度高， 计算时间稳定，没有收敛性问题等优点；尤其是当矩阵不变，而具有多个右端项或右端项不断变化 时，效率比迭代法更高。但是直接法程序实现的难度非常的大，其中所用到的算法通常深度抽象， 逻辑交错庞杂，即使是作者本人，在事隔几年后读自己的代码也废了不少精力，所以希望读者若对 这一领域感兴趣一定要有足够的耐心和毅力；若工作中并不涉及这一领域则建议跳过本章以免浪费 太多精力。稀疏直接求解方法繁多纷杂，但整体步骤大体一致： 1 顶点重排 2 符号分解 3 数值分解 4 三角回代 5 迭代改善 有时也将第1步和第2步合并在一起统称为符号分析，最后一步的迭代改善则是根据具体需求 决定是否使用，不是必须的。本章以求解大规模对称正定矩阵为例讲解开发快速求解大规模稀疏矩 阵的主要算法和相关优化技术。对称正定矩阵的求解技术已经发展的非常成熟，通常使用cholesky三 角分解方法，该方法数值稳定性好，且分解过程中无需选主元。在开始之前我们对一些符号进行约 定： adj(v) deg(v) : 顶点v的邻接节点的集合。 : 顶点v的度，也就是顶点为邻接节点的个数。reach(v): 顶点v的可达集，即通过顶点v可以到达的顶点的集合。snode( v0, v1, v2, …, vn )：表示组成超节点的所有顶点编号的集合，简写为snode。 innz(v) : 顶点v对应的波前中的非零元的行号或列号。3.1 基于quotient graph的符号分析 在进行真正的数值三角分解之前，我们需要事先确定存储分解因子所需要的内存大小，以避 免在分解的过程中动态内存分配。这一过程称之为符号分解或符号消元，这不仅可以提高程序的运 行效率，还可以节省很多内存并避免更多的内存碎片。假设矩阵A具有n个非零元，分解后的分解因 子L具有m个非零元(因此有m-n个非零元填入，通常mn&&n)，那么按照定义进行符号分解需要O(m)的存储空间。但是[1]中证明了，只需要O(n)大小的内存即 可完成符号分解，且这一大小是可以事先确定的。实现这一方法的模型就是quotient-graph。quotientgraph的执行流图如下：3.1.1 顶点重排 我们知道，三角分解过程中会产生很多填入元，如果直接进行分解，其填入元数量一般非常 多。当矩阵达到一定规模时，对存储空间的需求会大幅度暴涨，求解就会非常耗时和困难甚至无法 完成。所以我们需要事先对原矩阵的稀疏结构图中的顶点进行重排序处理。重排序的过程并不改变 原始图的拓扑结构，因而具有等价关系，下面两幅图展示了重排序的作用，左边是排序前的原始矩 阵，右边是排序后的矩阵：fig3.1fig3.2可以看到，如果直接进行分解，矩阵的所有元素都编程了非零元素；而重排序后进行分解非 零元素的数量不变，这虽然只是最理想的情况，但实际上通过顶点的重排序仍然能大幅度减少分解 过程中新产生的非零元素。稀疏矩阵顶点图的重排序算法通常都是NP问题，从而只能采用启发式算 法。目前最常用的算法有以缩减外形(或带宽)为目的算法，的如RCM算法； 以减少填入元为目的的类最小度算法(minimum degree)和用于分而治之的嵌套剖分(nested dissection)算法。虽然更窄的带宽通常意味着更少的填入元，但是其效果一般不如类最小度算法。嵌套 剖分目前比较成熟的技术是通过multilevel方法(类似多重网格方法中的思想)找到极小顶点分割子将整个 图剖分成多个独立的子图，然后再采用类最小度算法对各个子图分别进行局部优化处理(fig3.3)。类最小度算法常用的有最小度算法(MD),多重最小度算法(MMD)和近似最小度算法(AMD)。多数情况下性能最 高效果最好的是近似最小度算法，本书中我们仅采用最小度算法。 最消度算法的思想是当一个顶点的度比较小时，当消去该顶点时引入新的边的概率会更小或 是引入的新边的数量会更少，因此所有最小度算法都是从一个具有最小度的顶点开始。第一个具有 最小度的顶点可以随机的选择，也可以选择所有具有相同最小度的顶点中具有最大伪直径的顶点。 为了简单，我们遍历当前未消元的顶点并选择第一个顶点度最小的顶点。 最小度算法 pesudo-code do{ v=one vertex of graph with minimum degree elimination: delete all edges adjacency to ‘v’ add new edges between the adjacency verties of ‘v’(if non-existent edges between them) update graph: remove ‘v’ from graph }while(graph.num_verties&0) 最小度算法的实现有很多种方法，这里作者仅给出自己的实现（基于quotient graph模型）。以下代码都经过数百次针对不同布局和大小的矩阵的测试；但是按照严格的标准来说 ，数百次的测试远远不够，因此作者不保证代码中没有BUG，请读者慎用。 首先，需要选择一个具有最小度的顶点,这里我们称之为主元顶点： int search_pivot( const int* deg, const char* flag, int n ) { int i, d, p, for( i=0; i&n; ++i ){ if( flag[i]==0 ) } d=deg[i]; p=i; while( (++i)&n ){ if( flag[i]==0 ){ q=deg[i]; if( q&d ){ d=q; p=i; } } } } 然后需要得到这个顶点的可达集: int gather_reaches( aux_param_t* params, const int64* ptr, const int* nbor, int e ) { int64 p, int i, id, nen, n=0; p=ptr[e]; q=p+params-&len[e]; params-&flag[e]=1; while( p&q ){ id=nbor[p++]; if( params-&flag[id]==0 ){ if( params-&mark[id]==0 ){ params-&reach[n++]= params-&mark[id]=1; } } else { nen=gather_enbors( params, ptr, nbor, id ); for( i=0; i& ++i ){ id=params-&nbn[i];if( params-&mark[id]==0 ){ params-&reach[n++]= } params-&mark[id]=2; } } } } gather_reaches函数通过遍历主元顶点所有未消去的邻接节点以及所有以消去的邻接节点的邻接节点来 获得主元顶点的可达集;gather_enbors用来获取某个以消去顶点的邻接顶点： int gather_enbors( aux_param_t* params, const int64* ptr, const int* nbor, int e ) { int k, id, v, n=0; int64 p, if( params-&esize[e]&0 ) { v=e; for( k=0; k&params-&esize[e]+1; ++k ){ p=ptr[v]; q=p+params-&len[v]; while( p&q ){ id=nbor[p++]; if( params-&flag[id]==0 ){ params-&nbn[n++]= } } v=params-&elink[v]; } } else { p=ptr[e]; q=p+params-&len[e]; while( p&q ){ id=nbor[p++]; if( params-&flag[id]==0 ){ params-&nbn[n++]= } } } } 在消元过程中，可能有多个具有相同邻接顶点的以消去顶点，为了避免重复计算，需要对顶点进 行”吸收”,亦即将多个具有相同邻接节点的以消去顶点合并到其中一个超顶点。顶点吸收可以减少冗 余计算，同时也是基于quotientgraph模型计算时所必须的，否则内存空间将无处容纳顶点在消元过程中增加的边: void absorb_elements( int* nbor, const int64* ptr, aux_param_t* params, int nr, int e ) { int i, id, enb, n, n_deads, n_slots, sign, int64 ep, p=ptr[e], q=p+params-&len[e]; enb=0; for( enb=0; p&q; ++p ){ if( params-&flag[( id=nbor[p] )]!=0 ){ params-&nbn[enb++]= } } if( enb==0 ){ for( id=0; id& ++id ){ params-&mark[params-&reach[id]]=0; } } for( p=0, i=0; i& ++i ){ id=params-&nbn[i]; params-&spot[p++]= n=params-&esize[id]; params-&flag[id]=2; for( q=0; q&n; ++q ){ params-&spot[p++]=( id=params-&elink[id] ); }} v=e; params-&esize[e]=(int)p; for( i=0; i&p; ++i ){ id=params-&spot[i]; params-&elink[v]= v= } i=n_slots=sign=0; v=e; for(;;) { ++n_ ep=ptr[v]; p= q=ptr[v+1]; while( p&q ){ if( i&=nr ){ sign=1; } nbor[p++]=params-&reach[i++]; } params-&len[v]=(int)(p-ep); if( sign!=0 ) v=params-&elink[v]; } params-&esize[e]=--n_ for( params-&flag[e]=2, v=0; v& ++v ) { id=params-&reach[v]; if( params-&mark[id]!=2 ){ params-&mark[id]=0; } params-&mark[id]=0; p=ptr[id]; q=p+params-&len[id]; while( p&q ){ if( params-&flag[nbor[p]]==2 ){ nbor[p]=e; } ++p; } ep=p; i=0; while( (++p)&q ){ if( params-&flag[nbor[p]]!=2 ){ params-&spot[i++]=nbor[p]; }; } n_deads=(int)(q-(++ep)-i); if( n_deads&0 ){ for( p=0; p&i; ++p, ++ep ){ nbor[ep]=params-&spot[p]; } params-&len[id]-=n_ } } params-&flag[e]=1; } 最后我们需要更新顶点的度，首先我们需要计算出当前消元顶点的可达集（是对可达集中的顶点的 度进行更新），为此遍历消元顶点的邻接节点，如果其邻接节点是以消去的顶点，则遍历该顶点的 邻接顶点，如果其邻接顶点未被标记，则加入可达集；如果消元顶点的邻接顶点是未消去顶点，则 将其加入可达集（如果未标记的话）。顶点v的可达集可表示为 Reach(v) = adj(s)∪∑adj(e) 接着我们遍历可达集中的顶点，对于可达集中的每个顶点，遍历其邻接顶点，如果其邻接顶点为以 消去顶点，则将该顶点的顶点度加1，否则不变。如果其邻接顶点为非消去顶点，则将其顶点度加1, 如此进行，直至可达集中的所有顶点处理完成，下面是具体的实现代码： void update_degrees( aux_param_t* params, const int64* ptr, const int* nbor, int nr ) { int64 p, int k, d, s, n, i, for( k=0; k& ++k ) { s=params-&reach[k]; if( params-&deg[s]==0 ) d=0; p=ptr[s]; q=p+params-&len[s]; params-&mark[s]=1; do{id=nbor[p]; if( params-&flag[id]==0 ){ if( params-&mark[id]==0 ){ params-&spot[d++]= params-&mark[id]=1; } } else { n=gather_enbors( params, ptr, nbor, id ); for( i=0; i&n; ++i ){ id=params-&nbn[i]; if( params-&mark[id]==0 ){ params-&spot[d++]= params-&mark[id]=1; } } } }while( (++p)&q ); params-&deg[s]=d; params-&mark[s]=0; for( i=0; i&d; ++i ){ params-&mark[params-&spot[i]]=0; } } } 现在可以组装成完整的基于quotient-graph的最小度重排序函数了： void mdo( int* iperm, char* aux, int64* Ap, int* Ai, int n ) { aux_param_ int64 start, int d, i, e, params.flag =(char*) params.mark =params.flag+n; params.deg =(int*)( params.mark+n ); params.len =params.deg+n; params.esize =params.len+n; params.elink =params.esize+n; params.reach =params.elink+n; params.nbn =params.reach+n; params.spot =params.nbn+n; for( start=0, i=0; i&n; ++i ) { params.flag [i]=0; params.mark[i]=0; params.esize[i]=0; end=Ap[i+1]; d=(int)(end-start); params.deg[i]=d; params.len[i]=d; start= } for( i=0; i&n-1; ++i ) { e=search_pivot( params.deg, params.flag, n ); iperm[i]=e; m=gather_reaches( &params, Ap, Ai, e ); if( m==0 ) absorb_elements( Ai, Ap, &params, m, e ); update_degrees( &params, Ap, Ai, m ); } for( i=0; i&n; ++i ){ if( params.flag[i]==0 ) }iperm[n-1]=i; } aux是一个辅助内存数组，之所以没有在内部分配，是因为在求解大规模矩阵时，通常会对多 个子矩阵的顶点图进行重排序，这时候可以分配与参与计算的线程数量对等数量的aux数组，每个线 程中的aux数组的大小是该线程中所要计算的图中尺寸最大的图所对应的大小，这样aux数组便可被各 个线程中的多个顶点图复用，而不用频繁的开辟和释放操作。通过重排序得到了新的顶点编号后， 就可以建立消去树了。需要说明的是，在我们的最小度算法中使用的是外度更新，而非标准的顶点 度更新。顶点v的外度定义为由于考虑了额外更深层次顶点的影响，因此使用顶点外度的效果通常都要好于使用标准顶点度。 最小度算法有很多改进措施，第一个改进方法是将多个无区别顶点同时消去，这一方法也就 是MMD算法。所谓无区别顶点就是当两个不相邻顶点的邻接顶点集合与自身的集合相同的顶点，定 义为 u∪adj(u) &==& v∪adj(v), 用可达集可表示为 reach(v, S)∪v &==& reach(u,S)∪u 第二个改进方法是进行不完全的度更新，也就是不需要每次都对完整的邻接顶点集合的所有 顶点进行更新，而只更新满足一定条件的具有最小度的顶点，这一方法称为近似最小度(AMD)算法，A MD算法多数情况下的排序效果比MD和MMD算法好，同时算法的时间复杂度也大幅降低。 3.1.2 构建消去树 void build_etree( char* aux, int* parents, const int64* Cp, const int* Ci, int n ) { int64 a, b, p, q, *Rp, * int *Ri, *ancestor, i, k, Rp=(int64*) ptr=Rp+n+1; Ri=(int*)(ptr+n); ancestors=Ri+Cp[n]; memset( Rp, 0, n*sizeof(int64) ); p=Cp[0]; for( k=0; k&n; ++k ){ q=Cp[k+1]; while( p&q ){ ++Rp[Ci[p++]]; } } a=Rp[0]; Rp[0]=0; for( k=0; k&n; ++k ){ b=Rp[k+1]; Rp[k+1]=a; a+=b; } memcpy( ptr, Rp, n*sizeof(int64) ); p=Cp[0]; for(k=0; k&n; ++k ){ q=Cp[k+1];while( p&q ){ Ri[ptr[Ci[p++]]++]=k; } } for( k=0; k&n; ++k ) { parents[k]=-1; ancestors[k]=-1; for( p=Rp[k]; p&Rp[k+1]; ++p ) { i=Ri[p]; while((i&=0)&(i&k)){ next=ancestors[i]; ancestors[i]=k; if( next&0 ){ parents[i]=k; } i= } } } } 其中parents存储了消元过程中每个顶点之间的依赖关系，这个信息将在后面符号分解和寻找超 节点的过程中被用到。 3.1.4 符号分解 我们队矩阵A的稀疏结构X进行重新编码得到具有新的稀疏结构Y的矩阵B，Y与X中顶点之间的 拓扑关系并未改变，但是却具有更少的填入元。然后通过稀疏结构Y进行符号分解。符号分解的目的 是获取分解过程中分解因子中各个非零元素的位置信息以及整个分解因子所需要的存储空间大小， 从而避免进行实际的数值分解的过程中动态的计算索引以及动态的内存分配，这样将具有更高的效 率并避免内存碎片的积累。符号分解之所以有效是因为很多顶点对应的波前(后续章节会讲到)可能进 行多次更新，若果直接进行数值分解，就需要多次重复的计算索引以及很多低效的内存开辟释放操 作，严重影响效率。符号分解的代码同样基于quotient-graph,前面的一些函数会得到复用。 首先我们需要先计算得到分解因子中每列的非零元素数量，并最终得到整个分解因子所需要 的内存大小： int64 symbolic_eliminate( char* aux, int* Ln, const int64* Yp, const int* Yi, int n ) { int64 nA=Hp[n]&&1; int* Ai=(int*)(aux+ELIMINATION_BUFFER_BYTES(n)); int64* counter=(int64*)(Ai+nA); int64* Ap=counter+n+1; int64 p, q, c, int k, for( p=0, k=0; k&n; ++k ){ q=Yp[k+1]; counter[k]=q-p; p=q; } for( p=0, k=1; k&=n; ++k ){ q=Yp[k]; while( p&q ){ ++counter[Yi[p++]]; } } for( Ap[0]=0, k=0; k&n; ++k ){ Ap[k+1]=Ap[k]+counter[k]; counter[k]=Ap[k];} for( p=0, k=0; k&n; ++k ){ q=Yp[k+1]; while( p&q ){ Ai[counter[Yi[p++]]++]=k; } } for( p=0, k=0; k&n; ++k ){ q=Yp[k+1]; c=counter[k]; while( p&q ){ Ai[c++]=Yi[p++]; } } aux_param_ params.flag =(char*) params.mark =params.flag+n; params.len =(int*)( params.mark+n ); params.esize =params.len+n; params.elink =params.esize+n; params.reach =params.elink+n; params.nbn =params.reach+n; params.spot =params.nbn+n; for( p=0, k=0; k&n; ++k ){ params.flag [k]=0; params.mark [k]=0; params.esize[k]=0; q=Ap[k+1]; params.len[k]=(int)(q-p); p=q; } nnz=0; for( k=0; k&n; ++k ){ m=gather_reaches( &params, Ap, Ai, k ); if( m==0 ) Ln[k]=m; nnz+=m; absorb_elements( Ai, Ap, &params, n, k ); } } symbolic_eliminate的返回值即是分解因子非零元的数量，而分解因子每列非零元的数量保存在Ln 中,然后我们可以在分解前一次性的分配索引数组的内存，并利用Ln中的信息计算每列非零元素的行 号以及存储位置： void symbolic_decompose( char* aux, int* Li, const int64* Lp, const int* Sp, const int* nodes, const int* superIDs, const int64* Ap, const int* Ai, int n, int n_snodes ) { int64 ii, ii0, ii1, pp, qq, * char* int *spot, *a; int i k, v, d, sid, mark= pos=(int64*)(aux+n); spot=(int*)(pos+n_snodes); memset( mark, 0, n ); for( k=0; k&n_ ++k ){ v=nodes[Sp[k+1]-1]; ii0=Ap[v]; ii1=Ap[v+1]; pp=Lp[k]; while(ii0&ii1){ i=Ai[ii0++]; if(i&v){ Li[pp++]=Ai[ii0++]; } } pos[k]= } for( k=0; k&n_ ++k ) { ii0=Lp[k]; ii1=Lp[k+1]; d=(int)(ii1-ii0); if((k&0)&(d&1)){ a=&Li[ii1]; sort( a, a+d ); } for( ii=ii0; ii&ii1; ++ii ) { v=Li[ii]; sid=superIDs[v]; if( sid==k ) pp=Lp[sid]+1; qq=pos[sid]; c=0; while( pp&qq ){ v=Li[pp++]; mark[v]=1; spot[c++]=v; } for( pp=ii+1; pp&ii1; ++pp ){ v=Li[pp]; if( mark[v]==0 ){ Li[qq++]=v; } } for( i=0; i&c; ++i ){ mark[spot[i]]=0; } pos[sid]= } } } 这里我们是基于超节点进行的符号分解，因此在执行symbolic_decompose之前还需要构建超节点， 超节点方法的相关内容将在2.3节进行介绍。注意红色部分的代码，由于符号分解的过程中我们只记 录了正在进行分解的列的行号，但是为了以正确的顺序进行后续的分解，需要在此之前对其中的行 号（或列号）从小到大排序。现在整个符号分析已经完成，将其过程总结如下： 1 对顶点进行重排序。 2 利用重排序得到的新顶点编码构建消去树，得到消元过程中结点依赖关系的数组parents。3 利用新的顶点编码和和第2步中得到的parents信息寻找并构造超节点，并根据超节点内顶点编 码的连续性把超节点标记为静态或动态的。 4 利用新的顶点编码和第3步中得到的超节点信息进行符号分解确定整个分解过程所需要的内存 大小和非零元素的索引信息。3.2 多波前法 多波前(multifrontal)法发展于有限元中的波前法，其目的是提高稀疏矩阵分解过程中的并行性， 是目前求解大规模及超大规模稀疏矩阵普遍使用的高效且成熟的方法。下面我们以易于理解的方式 通过多个图的展示来说明多波前法的思想。 图4 第一步：波前{0,1,2,4}可以同时进行分解，波前0的结果对{2,6,7}进行更新， 波前1的结果对波前3进行更新，波前2和波前4的结果对波前5进行更新。 第二步：波前{3,5}可以同时进行分解，波前3的结果更新波前6。 第三步：对波前6进行分解，并用分解结果更新波前7。 第四步：对波前7进行分解。3.3 超结点法 超结点(supernodal)法也是用来提高计算效率的一种方法，但不同于多波前法同时对多个独立的 波前进行并行计算， 而是通过提高计算密度和缓存友好来提升效率。假设一个连续的顶点序列{ vi, vi+1, vi+2, …, vi+n, … }，如果 E1 parent(vi)=vi+1, parent(vi+1)=vi+2,…； 且 E2 num_adj(vi)+1=num_adg(vi+1), num_adj(vi+1)+1=num_adj(vi+2), … 那么将该顶点序列称之为一个超结点，如图所示： 图5 图中有5个超节点，分别为{0,1}，{2}，{3}，{4}，{5,6,7}。对超节点进行分解可以利用高效的稠密 矩阵计算提高计算效率。一般来说，获得超节点越少，表示消去树的高度越矮，计算密度越大，从 而效率也越高，因此可以通过一些方法来减少超节点的数量。第一种方法是超节点融合，亦即把超 节点当做普通的顶点，那么又可以对超节点进一步融合成更大的超节点；但是超节点融合会占用更 多的内存。第二种方法是在建立超节点时有限度的容忍，超节点中顶点的邻接节点的数量无需严格 满足E2，只要相差小于阀值即可加入超节点。第三种方法是扩展超节点的概念，超节点中的顶点无 需满足连续性，但E1条件仍然要满足，可将其称之为动态超节点。现实中更有效的方法是：将超节点分为静态超节点和动态超节点，同时配合第二种方法，这样不仅可以避免不必要的内存移动(因为 动态超节点在计算之前必须要将其移入一块连续的内存中)，还可以最大限度的提高计算的密度，且 内存的整体需求也一般会比较低。将超节点方法引入符号分解中还可以大幅降低存储索引所需的内 存以及减少符号分解所需的计算量，如图： 图6 图中包括两个动态超节点{0,2},{1,3}和一个静态超节点{4,5,6,7}。 容易看出，对于某个超节点snode{v0,v1,v2,…,vn},只需要存储其具有最小波前的顶点vn对应的那列 的非零元的行号，而其它顶点对应的波前中非零元的行号可以表示为: innz(v0) = { v1, v2, …, vn } + innz(vn) innz(v1) = { v2, v3, …, vn } + innz(vn) innz(v2) = { v3, v4, …, vn } + innz(vn) …… 从而整个超节点的非零元的行号（或列号）信息可以通过超节点中的包含的顶点的标号以及 最后一个顶点的非零元行号（或列号）完整的表示，亦即snode{ v0, v1, v2, …, vn } + innz(vn)，容易看出这样将使索引存储具有最小的空间需求。 建立超节点可以从正向寻找 (从第一列开始向后寻找)，也可以反向寻找(从最后一列开始向前寻找)，不过反向方法通常具有更好 的效果，这里分别给出