本文详细探讨了如何利用GPU进行高效编程，重点在于优化稀疏直接求解器的实现。通过一系列步骤和示例图片，展示了在GPU上处理大规模计算任务的策略和技巧，旨在提升计算速度和性能。
第三章 基于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 多波前+超节点方法的并行分解算法
小结
参考资料
第三章 基于GPU的稀疏直接求解器
前言
本章可能是所有章节中最难得了，如果节奏过快，对于即使有一定相关开发经验的人来说要彻底理解其中的内容也会相当困难，所以作者会给出相比其它章节更多的说明。第一章中类似的技术将会在本章得到应用，通过阅读本章，读者不仅可以进一步巩固第一章中学到的GPU计算的优化技术，还能了解到一些有趣的算法。在这个过程中读者还会更深的体会到从算法的伪代码到实际的可执行代码之间的难度。
3.0 稀疏矩阵直接解法概述
很多科学计算问题中都涉及到求解线性方程组，而这些方程组不仅规模很大且通常是稀疏的。解决这类问题通常使用迭代解法或直接解法。迭代解法一般内存空间需求小，程序实现比较简单；但是迭代解法的收敛速度依赖于矩阵的性态，如果矩阵过于病态，即使采用复杂的预处理技术(预处理越复杂，也意味着求解效率越低，内存空间需求越大，从而迭代解法优势也会逐渐丧失)也很难收敛甚至不收敛。直接解法的缺点是内存空间的需求比迭代解法大很多，但是其具有解的精度高，计算时间稳定，没有收敛性问题等优点；尤其是当矩阵不变，而具有多个右端项或右端项不断变化时，效率比迭代法更高。但是直接法程序实现的难度非常的大，其中所用到的算法通常深度抽象，逻辑交错庞杂，即使是作者本人，在事隔几年后读自己的代码也废了不少精力，所以希望读者若对这一领域感兴趣一定要有足够的耐心和毅力；若工作中并不涉及这一领域则建议跳过本章以免浪费太多精力。稀疏直接求解方法繁多纷杂，但整体步骤大体一致：
1 顶点重排
2 符号分解
3 数值分解
4 三角回代
5 迭代改善
有时也将第1步和第2步合并在一起统称为符号分析，最后一步的迭代改善则是根据具体需求决定是否使用，不是必须的。本章以求解大规模对称正定矩阵为例讲解开发快速求解大规模稀疏矩阵的主要算法和相关优化技术。对称正定矩阵的求解技术已经发展的非常成熟，通常使用cholesky三角分解方法，该方法数值稳定性好，且分解过程中无需选主元。在开始之前我们对一些符号进行约定：
adj(v) :顶点v的邻接节点的集合。
deg(v) : 顶点v的度，也就是顶点为邻接节点的个数。
reach(v) :顶点v的可达集，即通过顶点v可以到达的顶点的集合。
snode( v0, v1, v2, …) ：表示组成超节点的所有顶点编号的集合，简写为
snode。
innz(v) :顶点v对应的波前中的非零元的行号或列号。
3.1 基于
quotientgraph的符号分析
在进行真正的数值三角分解之前，我们需要事先确定存储分解因子所需要的内存大小，以避免在分解的过程中动态内存分配。这一过程称之为符号分解或符号消元，这不仅可以提高程序的运行效率，还可以节省很多内存并避免更多的内存碎片。假设矩阵A具有n个非零元，分解后的分解因子L具有m个非零元(因此有m-n个非零元填入，通常m-n>>n)，那么按照定义进行符号分解需要O(m)的存储空间。但是可以证明，只需要O(n)大小的内存即可完成符号分解，且这一大小是可以事先确定的。实现这一方法的模型就是
quotientgraph，我们后续符号分析中的顶点重排(
reordering)，构建消去树以及寻找超节点都是基于
quotientgraph的思想在O(n)的复杂度下完成的。
3.1.1 顶点重排
我们知道，三角分解过程中会产生很多填入元，如果直接进行分解，其填入元数量一般非常多。当矩阵达到一定规模时，对存储空间的需求会大幅度暴涨，求解就会非常耗时和困难甚至无法完成。所以我们需要事先对原矩阵的稀疏结构图中的顶点进行重排序处理。重排序的过程并不改变原始图的拓扑结构，因而具有等价关系，下面两幅图展示了重排序的作用，左边是排序前的原始矩阵，右边是排序后的矩阵：
fig3.1 fig3.2
可以看到，如果直接进行分解，矩阵的所有元素都编程了非零元素；而重排序后进行分解非零元素的数量不变，这虽然只是最理想的情况，但实际上通过顶点的重排序仍然能大幅度减少分解过程中新产生的非零元素。稀疏矩阵顶点图的重排序算法通常都是NP问题，从而只能采用启发式算法。目前最常用的算法有以缩减外形(或带宽)为目的算法，的如
RCM算法；以减少填入元为目的的类最小度算法(
minimum degree)和用于分而治之的嵌套剖分(
nested dissection)算法。虽然更窄的带宽通常意味着更少的填入元，但是其效果一般不如类最小度算法。嵌套剖分目前比较成熟的技术是通过
multilevel方法(类似多重网格方法中的思想)找到极小顶点分割子将整个图剖分成多个独立的子图，然后再采用类最小度算法对各个子图分别进行局部优化处理(fig3.3)。类最小度算法常用的有最小度算法(
MD),多重最小度算法(
MMD)和近似最小度算法(
AMD)。多数情况下性能最高效果最好的是近似最小度算法，本书中我们仅采用最小度算法。、
（白色部分为0元素）
最消度算法的思想是当一个顶点的度比较小时，当消去该顶点时引入新的边的概率会更小或是引入的新边的数量会更少，因此所有最小度算法都是从一个具有最小度的顶点开始。第一个具有最小度的顶点可以随机的选择，也可以选择所有具有相同最小度的顶点中具有最大伪直径的顶点。为了简单，我们遍历当前未消元的顶点并选择第一个顶点度最小的顶点。
最小度算法
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 betweenthem)
update graph: remove ‘v’ from graph
}while(graph.num_verties>0)
最小度算法的实现有很多种方法，这里作者仅给出自己的实现（基于quotientgraph模型）。以下代码都经过数百次针对不同布局和大小的矩阵的测试；但是按照严格的标准来说，数百次的测试远远不够，因此作者不保证代码中没有BUG，请读者慎用。
首先，需要选择一个具有最小度的顶点,这里我们称之为主元顶点：
int
search_pivot(
const
int
* deg,
const
char
* flag,
int
n )
{
int
i, d, p, q;
for
( i=0; i<n; ++i ){
if
(flag==0 ) break; }
d=deg; p=i;
while( (++i)<n ){ if( flag==0 ){ q=deg; if( q<d){ d=q; p=i; } } }
return p;
}
然后需要得到这个顶点的可达集:
int gather_reaches( aux_param_t* params,const int64* ptr, const int* nbor, int e )
{
int64 p, q;
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++]=id;params->mark[id]=1;
}
} else {
nen=gather_enbors(params, ptr, nbor, id );
for( i=0; i<nen; ++i ){
id=params->nbn;
if( params->mark[id]==0 ){ params->reach[n++]=id; }
params->mark[id]=2;
}
}
}
return n;
}
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, q;
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++]=id; }
}
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++]=id; }
}
}
return n;
}
在消元过程中，可能有多个具有相同邻接顶点的以消去顶点，为了避免重复计算，需要对顶点进行”吸收”,亦即将多个具有相同邻接节点的以消去顶点合并到其中一个超顶点。顶点吸收可以减少冗余计算，同时也是基于quotient-graph模型计算时所必须的，否则内存空间将无处容纳顶点在消元过程中增加的边:
void absorb_elements( <
最低0.47元/天 解锁文章
949
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
