第三章 基于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)大小的内存即可完成符号分解，且这一大小是可以事先确定的。实现这一方法的模型就是[b][i]quotientgraph[/i][/b]，我们后续符号分析中的顶点重排([b][i]reordering[/i][/b])，构建消去树以及寻找超节点都是基于[b][i]quotientgraph[/i][/b]的思想在O(n)的复杂度下完成的。
3.1.1 顶点重排
我们知道，三角分解过程中会产生很多填入元，如果直接进行分解，其填入元数量一般非常多。当矩阵达到一定规模时，对存储空间的需求会大幅度暴涨，求解就会非常耗时和困难甚至无法完成。所以我们需要事先对原矩阵的稀疏结构图中的顶点进行重排序处理。重排序的过程并不改变原始图的拓扑结构，因而具有等价关系，下面两幅图展示了重排序的作用，左边是排序前的原始矩阵，右边是排序后的矩阵：
[attach]3511[/attach] [attach]3512[/attach]
fig3.1 fig3.2
可以看到，如果直接进行分解，矩阵的所有元素都编程了非零元素；而重排序后进行分解非零元素的数量不变，这虽然只是最理想的情况，但实际上通过顶点的重排序仍然能大幅度减少分解过程中新产生的非零元素。稀疏矩阵顶点图的重排序算法通常都是NP问题，从而只能采用启发式算法。目前最常用的算法有以缩减外形(或带宽)为目的算法，的如[b][i]RCM[/i][/b]算法；以减少填入元为目的的类最小度算法([b][i]minimum degree[/i][/b])和用于分而治之的嵌套剖分([b][i]nested dissection[/i][/b])算法。虽然更窄的带宽通常意味着更少的填入元，但是其效果一般不如类最小度算法。嵌套剖分目前比较成熟的技术是通过[b][i]multilevel[/i][/b]方法(类似多重网格方法中的思想)找到极小顶点分割子将整个图剖分成多个独立的子图，然后再采用类最小度算法对各个子图分别进行局部优化处理(fig3.3)。类最小度算法常用的有最小度算法([b][i]MD[/i][/b]),多重最小度算法([b][i]MMD[/i][/b])和近似最小度算法([b][i]AMD[/i][/b])。多数情况下性能最高效果最好的是近似最小度算法，本书中我们仅采用最小度算法。、
[attach]3513[/attach]
（白色部分为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 ){ ifi [/i]break; }
[i] d=deg; p=i;
whilei{ [/i]ifi{ q=deg; ifi{ d=q; p=i; } } }[/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;
whilei{[/i]
id=nbor[p++];
ifi{[/i]
ifi{ [/i]
params->reach[n++]=id;params->mark[id]=1;
}
} else {
nen=gather_enbors(params, ptr, nbor, id );
for( i=0; i<nen; ++i ){
[i] id=params->nbn;
ifi{ params->reach[n++]=id; }[/i]
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;
ifi[/i]
{
v=e;
for( k=0; kesize[e]+1;++k ){
p=ptr[v];q=p+params->len[v];
whilei{ [/i]
id=nbor[p++];
ifi{ params->nbn[n++]=id; }[/i]
}
v=params->elink[v];
}
}
else
{
p=ptr[e];q=p+params->len[e];
whilei{ [/i]
id=nbor[p++];
ifi{ params->nbn[n++]=id; }[/i]
}
}
return n;
}
在消元过程中，可能有多个具有相同邻接顶点的以消去顶点，为了避免重复计算，需要对顶点进行”吸收”,亦即将多个具有相同邻接节点的以消去顶点合并到其中一个超顶点。顶点吸收可以减少冗余计算，同时也是基于quotient-graph模型计算时所必须的，否则内存空间将无处容纳顶点在消元过程中增加的边:
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, v;
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++]=id; }
}
ifi{ [/i]
for( id=0; id<nr; ++id ){ params->mark[params->reach[id]]=0; } return;
}
for( p=0, i=0; i<enb; ++i ){
[i] id=params->nbn;
params->spot[p++]=id;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]( i=0; i<p; ++i ){ id=params->spot; params->elink[v]=id; v=id;}
i=n_slots=sign=0; v=e;
fori[/i]
{
++n_slots;ep=ptr[v]; p=ep; q=ptr[v+1];
whilei{[/i]
ifi{ sign=1; [/i]break; }
[i] nbor[p++]=params->reach;
}
params->len[v]=(int)(p-ep);
ifi [/i]break;
v=params->elink[v];
}
params->esize[e]=–n_slots;
for( params->flag[e]=2, v=0; v<nr; ++v )
{
id=params->reach[v];
ifi{ params->mark[id]=0; [/i]continue; }
params->mark[id]=0;p=ptr[id]; q=p+params->len[id];
whilei{ [/i]ifi{nbor[p]=e; [/i]break; } ++p; }
ep=p; i=0;
whilei{ [/i]
ifi{ [/i]
[i]params->spot=nbor[p];
}
}
n_deads=(int)(q-(++ep)-i);
ifi{ [/i]
for( p=0; p<i; ++p, ++ep ){ nbor[ep]=params->spot[p]; }
params->len[id]-=n_deads;
}
}
params->flag[e]=1;
}
最后我们需要更新顶点的度，首先我们需要计算出当前消元顶点的可达集（是对可达集中的顶点的度进行更新），为此遍历消元顶点的邻接节点，如果其邻接节点是已消去的顶点，则遍历该顶点的邻接顶点，如果其邻接顶点未被标记，则加入可达集；如果消元顶点的邻接顶点是未消去顶点，则将其加入可达集（如果未标记的话）。顶点v的可达集可表示为
[b][i]reach(v)= adj(s)∪∑adj(e)[/i][/b]
接着我们遍历可达集中的顶点，对于可达集中的每个顶点，遍历其邻接顶点，如果其邻接顶点为已消去顶点，则将该顶点的顶点度加1，否则不变。如果其邻接顶点为非消去顶点，则将其顶点度加1,如此进行，直至可达集中的所有顶点处理完成，下面是具体的实现代码：
voidupdate_degrees( aux_param_t* params, const int64* ptr, const int* nbor, int nr )
{
int64 p, q;
int k, d, s, n, i, id;
for( k=0; k<nr; ++k )
{
s=params->reach[k];
ifi [/i]continue;
d=0;p=ptr[s]; q=p+params->len[s]; params->mark[s]=1;
do{
id=nbor[p];
ifi{[/i]
ifi{ [/i]
params->spot[d++]=id;
params->mark[id]=1;
}
}else {
n=gather_enbors(params, ptr, nbor, id );
for( i=0; i<n; ++i ){
[i] id=params->nbn;
ifi{ [/i]
params->spot[d++]=id;params->mark[id]=1;
}
}
}
}whilei;[/i]
params->deg[s]=d;params->mark[s]=0;
for[i]( i=0; i<d; ++i ){ params->mark[params->spot]=0; }
}
}
现在可以组装成完整的基于quotient-graph的最小度重排序函数了:
void mdo( int* iperm, char* aux, int64* Ap, int* Ai, int n )
{
aux_param_t params;
int64 start, end;
int d, i, e, m;
params.flag =(char*)aux;
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 )
{
[i] params.flag=0;
[i] params.mark=0;
[i] params.esize=0;
[i] end=Ap;
d=(int)(end-start);
[i] params.deg=d;
[i] params.len=d;
start=end;
}
for( i=0; i<n-1; ++i )
{
e=search_pivot(params.deg, params.flag, n );
[i] iperm=e;
m=gather_reaches(¶ms, Ap, Ai, e );
ifi [/i]continue;
absorb_elements(Ai, Ap, ¶ms, m, e );
update_degrees(¶ms, Ap, Ai, m );
}
for( i=0; i<n; ++i ){ ifi [/i]break; }
iperm[n-1]=i;
}
aux是一个辅助内存数组，之所以没有在内部分配，是因为在求解大规模矩阵时，通常会对多个子矩阵的顶点图进行重排序，这时候可以分配与参与计算的线程数量对等数量的aux数组，每个线程中的aux数组的大小是该线程中所要计算的图中尺寸最大的图所对应的大小，这样aux数组便可被各个线程中的多个顶点图复用，而不用频繁的开辟和释放操作。通过重排序得到了新的顶点编号后，就可以建立消去树了。需要说明的是，在我们的最小度算法中使用的是外度更新，而非标准的顶点度更新。由于考虑了额外更深层次顶点的影响，因此使用顶点外度的效果通常都要好于使用标准顶点度。
最小度算法有很多改进措施，第一个改进方法是将多个无区别顶点同时消去，这一方法也就是MMD算法。所谓无区别顶点就是当两个不相邻顶点的邻接顶点集合与自身的集合相同的顶点，定义为
u∪adj(u) <==> v∪adj(v)
用可达集可表示为
reach(v, S)∪v <==> reach(u,S)∪u
第二个改进方法是进行不完全的度更新，也就是不需要每次都对完整的邻接顶点集合的所有顶点进行更新，而只更新满足一定条件的具有最小度的顶点，这一方法称为近似最小度(AMD)算法，AMD算法多数情况下的排序效果比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,*ptr;
int *Ri,*ancestor, i, k, next;
Rp=(int64*)aux;
ptr=Rp+n+1;
Ri=(int*)(ptr+n);
ancestors=Ri+Cp[n];
memset( Rp, 0, n*sizeofi );[/i]
p=Cp[0];
for( k=0;k<n; ++k ){
q=Cp[k+1];
whilei{ ++Rp[Ci[p++]]; }[/i]
}
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*sizeofi );[/i]
p=Cp[0];
for(k=0;k<n; ++k ){
q=Cp[k+1];
whilei{ Ri[ptr[Ci[p++]]++]=k; }[/i]
}
for( k=0;k<n; ++k )
{
parents[k]=-1;ancestors[k]=-1;
for( p=Rp[k]; p<Rp[k+1]; ++p )
{
i=Ri[p];
whilei{[/i]
[i]next=ancestors;
[i] ancestors=k;
ifi{ parents=k; }
i=next;
}
}
}
}
其中parents存储了消元过程中每个顶点之间的依赖关系，这个信息将在后面符号分解和寻找超节点的过程中被用到。
3.1.4 符号分解
我们队矩阵A的稀疏结构X进行重新编码得到具有新的稀疏结构Y的矩阵B，Y与X中顶点之间的拓扑关系并未改变，但是却具有更少的填入元。然后通过稀疏结构Y进行符号分解。符号分解的目的是获取分解过程中分解因子中各个非零元素的索引以及整个分解因子所需要的存储空间大小，从而避免在数值分解的过程中动态的计算索引以及动态的内存分配，这样将具有更高的效率并避免内存碎片的积累。符号分解之所以有效是因为很多顶点对应的波前(后续章节会讲到)可能进行多次更新，如果直接进行数值分解，就需要多次重复的计算索引以及很多低效的内存分配释放操作，严重影响效率。符号分解的代码同样基于[b][i]quotient-graph[/i][/b],前面的一些函数会得到复用。
首先我们需要先计算得到分解因子中每列非零元素的数量，并最终得到整个分解因子所需要的内存大小：
int64 symbolic_eliminate( char* aux, int* Ln, const int64* Yp, const int* Yi, int n )
{
int64nA=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, nnz;
int k, m;
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];
whilei{ ++counter[Yi[p++]]; }[/i]
}
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];
whilei{ Ai[counter[Yi[p++]]++]=k; }[/i]
}
for( p=0,k=0; k<n; ++k ){
q=Yp[k+1];
c=counter[k];
whilei{ Ai[c++]=Yi[p++]; }[/i]
}
aux_param_t params;
params.flag =(char*)aux;
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(¶ms, Ap, Ai, k );
ifi [/i]continue;
Ln[k]=m; nnz+=m;
absorb_elements( Ai, Ap,¶ms, n, k );
}
return nnz;
}
[b][i]symbolic_eliminate[/i][/b]的返回值是分解因子非零元的数量，而分解因子每列非零元的数量保存在[b][i]Ln[/i][/b]中,然后我们可以在分解前一次性的分配索引数组的内存，并利用[b][i]Ln[/i][/b]获取每列非零元素的行号以及存储位置：
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, intn_snodes )
{
int64 ii, ii0, ii1, pp,qq, *pos;
char*mark;
int *spot,*a;
int i k,v, d, sid, c;
mark=aux;
pos=(int64*)(aux+n);
spot=(int*)(pos+n_snodes);
memset( mark, 0, n );
for( k=0;k<n_snodes; ++k )
{
v=nodes[Sp[k+1]-1];
ii0=Ap[v];
ii1=Ap[v+1];
pp=Lp[k];
whilei{ [/i]
i=Ai[ii0++];
ifi{ Li[pp++]=Ai[ii0++]; }[/i]
}
pos[k]=pp;
}
for( k=0;k<n_snodes; ++k )
{
ii0=Lp[k];
ii1=Lp[k+1];
d=(int)(ii1-ii0);
ifi{[/i]
a=&Li[ii1];
sort( a,a+d );
}
for( ii=ii0; ii<ii1; ++ii )
{
v=Li[ii];
sid=superIDs[v];
ifi [/i]continue;
pp=Lp[sid]+1;
qq=pos[sid];
c=0;
whilei{[/i]
v=Li[pp++];mark[v]=1; spot[c++]=v;
}
for( pp=ii+1; pp<ii1; ++pp ){
v=Li[pp];
ifi{ Li[qq++]=v; }[/i]
}
for[i]( i=0; i<c; ++i ){ mark[spot]=0; }
pos[sid]=qq;
}
}
}
这里我们是基于超节点进行的符号分解，因此在执行symbolic_decompose之前还需要构建超节点，超节点方法的相关内容将在2.3节进行介绍。注意红色部分的代码，由于符号分解的过程中我们只记录了正在进行分解的列的行号，但是为了以正确的顺序进行后续的分解，需要在此之前对其中的行号（或列号）从小到大排序。现在整个符号分析已经完成，将其过程总结如下：
1 对顶点进行重排序。
2 利用重排序得到的新顶点编码构建消去树，得到消元过程中结点依赖关系的数组parents。
3 利用新的顶点编码和第2步中得到的parents信息寻找并构造超节点，并根据超节点内顶点编码的连续性把超节点标记为静态超节点或动态超节点的。
4 利用新的顶点编码和第3步中得到的超节点信息进行符号分解确定整个分解过程所需要的内存大小和非零元素的索引信息。
3.2 多波前法（multifrontalmethod）
多波前法发展于有限元中的波前法，其目的是提高稀疏矩阵分解过程中的并行性，是目前求解大规模及超大规模稀疏矩阵普遍使用的高效且成熟的方法。下面我们以易于理解的方式通过图的展示来说明多波前法的思想（假设已经得到了分解后的稀疏结构）。
[attach]3514[/attach]
[attach]3515[/attach]
第一步：波前[b][i]{0,1,2}[/i][/b]可以同时进行分解，波前0的结果对[b][i]{3，4,7}[/i][/b]进行更新，
波前1的结果对波前{4，6}进行更新，波前2对波前{5，6}进行更新。
第二步：对波前[b][i]{3}[/i][/b]进行分解，然后更新波前4。
第三步：对波前[b][i]{6，7}[/i][/b]进行并行分解。
3.3 超结点法（supernodalmethod）
超结点法也是用来提高计算效率的一种方法，但不同于多波前法同时对多个独立的波前进行并行计算， 而是通过提高计算密度和缓存友好来提升效率。假设一个连续的顶点序列{ 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), …
那么将该顶点序列称之为一个超结点，如图所示：
[b] [/b] [attach]3516[/attach]
图中有5个超节点（红色cell代表分解过程中生成的填入元素，其它代表原始的稀疏结构），分别为{1}， {2}, {0，3, 6,7}，{4，5}。如果按照超节点进行计算，其顺序为：
a. 首先对超节点{1},{2},{4,5}进行分解（由于没有依赖关系，所以三个超节点可同时进行分解）
b. 使用分解得到的波前矩阵对3，6，7列进行更新
c. 对超节点{0,3,6,7}进行分解
对超节点进行分解可以利用高效的稠密矩阵计算提高计算效率。一般来说，超节点越少，则消去树的高度越矮，计算密度越大，从而效率也越高，因此可以通过一些方法来减少超节点的数量。第一种方法是超节点融合，亦即把超节点当做普通的顶点，那么又可以对超节点进一步融合成更大的超节点；但是超节点融合会占用更多的内存。第二种方法是在建立超节点时有限度的容忍，超节点中顶点的邻接节点的数量无需严格满足E2，只要相差小于阀值即可加入超节点，仍以上图为例，可以将超节点{1}和{0,3,6,7}合并为一个超节点：
{0，1，3，6，7}
这样虽然浪费了一个元素的空间，但是却可以减少超节点的数量（意味着消去树层级更少，也因此减少了消去过程中的同步开销，同时计算密度更大从而具有更高的效率）。
第三种方法是扩展超节点的概念，超节点中的顶点无需满足连续性，但E1条件仍然要满足，可将其称之为动态超节点。现实中更有效的方法是：将超节点分为静态超节点和动态超节点，同时配合第二种方法，这样不仅可以避免不必要的内存移动(因为动态超节点在计算之前必须要将其移入一块连续的内存中)，还可以最大限度的提高计算的密度，且内存的整体需求也一般会比较低。将超节点方法引入符号分解中还可以大幅降低存储索引所需的内存并减少符号分解所需的计算量。
对于上图来说，{1}, {2}, {4,5} 是静态超节点，{0,3,6,7}是动态超节点。
容易看出，对于某个超节点[b][i]snode{v0,v1,v2,…,vn}[/i][/b],只需要存储其具有最小波前的顶点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)，容易看出这样将使索引存储具有最小的空间需求。
建立超节点可以从正向寻找 (从第一列开始向后寻找)，也可以反向寻找(从最后一列开始向前寻找)，不过反向方法通常具有更好的效果，这里分别给出两种实现，在使用中可以从前向和后向方法中选择结果最好的。
寻找超结点的反向方法：
int find_supernodals_backward( char* aux, int* Sp, int*verties, char* flags, const int64* Lp, const int* Li, int n )
{
int64 ii0, ii1, ii2, ii3;
[i]int[/i][i]* child=([/i][i]int[/i][i]*)aux;[/i]
int*id=child+n;
int i, k,n_snodes, snode_size, p, q, parent;
char b;
for( k=0;k<n; ++k ){ child[k]=-1; }
for( k=0;k<n-1; ++k )
{
ii0=Lp[k];
ii1=Lp[k+1];
parent=Li[ii0+1];
ii2=Lp[parent];
ii3=Lp[parent+1];
ifi{[/i]
child[parent]=k;
}
}
for( i=0,n_snodes=0, k=n-1; k>=0; --k )
{
ifi[/i]
{
snode_size=1;
[i] id=k;
p=child[k];
child[k]=-2;
whilei{[/i]
++snode_size;
[i] id=p;
q=child[p];
child[p]=-2;
p=q;
}
Sp[n_snodes++]=snode_size;;
}
}
for[i]( i=0;i<n; ++i ){ verties=id[n-i-1]; }
for[i]( i=0;i<n_snodes; ++i ){ child=Sp[n_snodes-i-1]; }
for[i](Sp[0]=0, i=0; i<n_snodes; ++i ){ Sp[i+1]=Sp[i]+child; }
for( k=0;k<n_snodes; ++k ){
char b=0;
for( i=Sp[k]; i<Sp[k+1]; ++i ){
ifi{ b=1; [/i]break; }
}
flags[k]=b;
}
returnn_snodes;
}
寻找超结点的正向方法：
int find_supernodals_forward( char* mark, int* Sp, int*verties, char* flags, const __int64* Lp, const int* Li, int n )
{
memset( mark, 0, n );
int64 ii0,ii1, ii2, ii3;
int i, k,parent, n_snodes, snode_size;
char b;
const int imax=n-1;
for(Sp[0]=0, i=0, n_snodes=0, k=0; k<n; ++k )
{
ifi [/i]continue;
snode_size=1;
mark[k]=1;
[i] verties=k;
ii0=Lp[k]; ii1=Lp[k+1];
parent=(k<imax)?Li[ii0+1]:n;
whilei{ [/i]
ii2=Lp[parent];ii3=Lp[parent+1];
ifi [/i]break;
++snode_size;
mark[parent]=1;
[i] verties=parent;
parent=(parent<imax)?Li[ii2+1]:n;
ii0=ii2;ii1=ii3;
}
Sp[n_snodes+1]=Sp[n_snodes]+snode_size;++n_snodes;
}
for( k=0;k<n_snodes; ++k ){
b=0;
for( i=Sp[k]; i<Sp[k+1]; ++i ){
ifi{ b=1; [/i]break; }
}
flags[k]=b;
}
returnn_snodes;
}
find_supernodal_*函数的返回值时超节点的数量；verties的大小是顶点的数量，每个超节点所包含的顶点编号在verties中是连续存储的；Sp中包含了超节点的第一个顶点在数组中的索引位置以及超节点的大小信息，第i个超节点中的第k个顶点为verties[Sp[i]+k], 其中k<Sp[i+1]。可以如下遍历第i个超节点中的顶点：
[i]p=Sp[i]；q=Sp;
do{
v=verties[p];
…
}while((++p)<q);
flags中则标记了超节点的静态和动态属性，值为0表示静态超节点，值为1表示动态超节点。
3.4 超结点的分解
对超结点的分解可以方便的利用高效的稠密矩阵算法，下图是超结点分解的示意图
由于大部分的计算都集中在syrk,因此我们重点对syrk的实现进行优化，第一章中的矩阵乘法所用到的技术可以得到复用，只不过矩阵B变成了矩阵A的转置。虽然可以通过直接使用*gemm来实现syrk操作，但是这种方法还是比较低效，因为每一步的计算量通常并不足够大，可能无法充分利用硬件的计算资源；而且syrk比gemm具有更少的全局内存访问量。另一个需要自己开发这些操作而不使用cublas等通用库的原因是它们的实现中syrk操作都是基于完整尺寸的矩阵，对于稀疏矩阵求解这种内存非常宝贵的计算来说过于浪费。因此我们实现自己的syrk操作通过对称存储且只需要大约一半的存储空间，其数据布局如下图所示：
[attach]3520[/attach]
3.5 基于多波前+超结点方法的并行分解
可以将多波前法和超节点方法结合使用以最大化并行度，每个超节点作为一个超波前(如无特别说明，我们将超节点波前统一简称为波前)，多个波前可以同时进行计算。在安排超节点的计算顺序时，除了根据已经建立的消去树外，还要考虑当前的计算资源，包括GPU的数量，每个GPU上得设备内存的大小以及主机内存的大小。我们需要在求解性能和存储空间之间取得一个最佳的平衡。因此每个设备上可以允许同时进行计算的波前需要设置一个上限以避免在消去树的某一层中的波前对存储空间的需求超出已有的内存资源。同时当有多个GPU设备时还需要综合考虑如何将计算和所需要的设备内存平均到每个设备上。通常对消去树的深度遍历可以最小化存储资源使用量，而宽度遍历则可最大化并行度。另一个需要考虑的就是当矩阵规模很大时，内存和显存往往只能容下小部分数据，这时候就需要考虑核外求解，也就是将已经参与完计算的波前数据写入磁盘。但是后续的计算可能会要求载入已经写入磁盘的波前数据，那么在确定波前的消去顺序时又需要考虑如何最小化磁盘读写操作。确定各个波前消去顺序的方法很多，但是波前消去的顺序的多选择性仅限于消去树同一层的波前之间，各个层之间的波前的消去仍需按顺序进行，因此可以很容易的根据自己的设想或需求设计相关算法，这里不再叙述。下面给出基于多波前+超节点方法的分解算法：
idev= current_cpu_tid : 当前设备的id
A: 待分解的矩阵
U[idev]: 第idev设备上得波前更新矩阵，用来存储每个波前的syrk的计算结果
for( level = 0; level<etree.n_levels; ++level )
{
start_snode =etree.snode_list[level][idev]
end_snode=etree.snode_list[level][idev+1]
get_front( front[idev],snode, A )
for( snode=start_snode+1;snode<end_snode; ++snode ){
if(next_snode.flag == dynamic ){
async{
assembly_front(front[idev].hFactor, next_snode, A )
copy(front[idev], next_snode )
}
} else {
front.[idev].hFactor=A[snode.Ap[next_snode.id]: snode.Ap[next_snode.id+1]];
async_copy(front, next_snode )
}
front[idev].dU= snode.factor( front[idev].dF )
locked_extend_add(A, front[idev], snode )
if((compute_mode==out_core)&&(snode.have_no_parents)&&(snode!=last)){
async_output_to_disk(snode )
}
}
}
上述算法是基于多GPU的多波前+超节点方法的实现，使用了静态+动态超节点分类技术。根据消去树的指导，逐层求解各个层内的波前(亦即超节点波前)。开始分解时整个待分解矩阵存储在内存或磁盘中，然后在每一层计算的开始处使用异步传输将下一个要计算的波前数据从主机内存到设备内存的复制过程与当前波前的计算重叠进行，从而可以有效隐藏CPU和GPU之间进行数据传输延迟，大部分时候这个延迟基本可以完全隐藏掉。然后通过判断当前超节点的标志，如果是静态超节点，只需直接进行原地计算；如果是动态超节点，则要先组装超节点再进行计算。计算的结果分为两部分，一个是波前本身的分解因子F，另一个是当前波前对其它波前的更新矩阵U。由于可能有多个波前对同一个波前的计算有贡献，因此我们需要将extend_add操作锁定以保证计算不发生写冲突。最后需要根据一定的条件决定是否将以完成计算的波前写入磁盘，如果设定的计算模式是核外求解模式并且没有后续的波前计算需要用到当前的波前，则将数据写入磁盘。同时还要注意的是第三个判别条件，如果该波前是整个消去树中的最后一个波前，那么考虑到后续需要进行回代求解，因此从效率上考虑，就无需再写回磁盘了。
本章小结
这一章我们介绍并给出了稀疏直接求解器的具体实现以及先关的优化技术，这其中的很多技术都可以直接或是稍加变通的用在很多数值算法中。特别是超节点方法特别适合用GPU进行计算，但是也需要通过各种优化尽量减少CPU和GPU之间的数据通信带来的开销，比如：
通过对超节点的静态和动态分类减少数据移动的开销，最大化原地分解。
通过动态并行对完整的三角型超节点进行分解以减少和主机通信的开销。
通过使用内存映射和点对点传输降低CPU和GPU之间以及GPU和GPU之间的数据复制开销。
通过对超节点波前的消去顺序进行调度排序使得每个波前进出内存的次数最小化。
使用基于双全局内存和双主机内存的双缓冲技术将待计算节点/波前的传输和当前节点/波前的计算重叠进行，从而充分利用CPU和GPU之间的异步并发模式。
参考资料
1《Direct Methods for Sparse LinearSystem》，Timothy A. Davis
2 An Approximate Minimum Degree Ordering Algorithm，Patrick R. Amestoy Timothy A.Davisy Iain S. Du_z SIAM J, 1996
3 The Computational Complexity of theMinimum Degree Algorithm，P. Heggernesy S. C. Eisenstatz G. Kumfertx A. Pothen，2000
4 [i]Theevolution of the minimum degree ordering algorithm[/i], SIAM Review，1989
5 J.A. GEORGE AND J. W. H. LIU, [i]A quotient graph model for symmetricfactorization[/i],
in SparseMatrix Proceedings 1978, I. S. Duff and G. W. Stewart, eds., SIAM
Publications,1978
6 A NESTED DISSECTION APPROACH TO SPARSE MATRIXPARTITIONING FOR PARALLEL COMPUTATIONS, ERIK G.BOMAN_ AND MICHAEL M. WOLF
7 Minimum Degree Reordering Algorithms: A Tutorial Stephen Ingram University of British Columbia