Transcript 大规模稀疏矩阵并行计算
大规模稀疏矩阵
并行计算
李修宇
QQ：295553381
4/8/2015
1
主流求解方法
• 直接法
o GAUSS消去法
o 波前法
o 多波前法
• 迭代法
o 经典迭代法
• Jacobi、SOR、SSOR
o 投影方法
• CG、GMRES
o 预处理技术
• 不完全分解预处理条件
o 代数多重网格技术
大规模稀疏矩阵并行计算
4/8/2015
2
矩阵性质对求解的影响
性质
• 非零元的分布
o 带状分布
o 按块分布
o ……
影响
• 矩阵的存储方式
• 求解方法的选择
• 正定性
• 求解速度
• 对称性
• ……
大规模稀疏矩阵并行计算
4/8/2015
3
直接法
• 矩阵图重排：一般分为两大类，带宽缩减算法（也常称为
外形缩减）和区域分解算法，应用较多的带宽缩减算法
CM，RCM，GPS，Rosen算法。一般建议多重方法结合
使用：全局方法的全局平衡性、局部方法的局部最优特性。
• 符号分解：确定非零元结构以及相应的消元索引，以便在
实际数值分解前确定所需存储资源大小，避免数值分解中
动态分配存储空间和复杂的索引策略。
• 构建消去树(elimination tree)：确定分解节点之间的分解
依赖，即确定分解的顺序并构成并行分解的层次结构。
大规模稀疏矩阵并行计算
4/8/2015
4
直接法
• 数值分解：利用符号分解得到的非零元结构和索引沿消去
树路径进行分解。
• 回代求解：包括前向（forward）和后向（backward）
回代，可先构建消去依赖树或顶点着色技术实现并行回代
求解。
• 在有限元领域应用最广的直接求解方法常使用带宽缩减或
多区域分解的多波前法（multifrontal）。
大规模稀疏矩阵并行计算
4/8/2015
5
对称正定矩阵的求解
大规模稀疏矩阵并行计算
4/8/2015
6
对称矩阵的不完全分解
大规模稀疏矩阵并行计算
4/8/2015
7
代数多重网格法
• V-Cycle AMG（V循环多重网格法）
• W-Cycle AMG（W循环多重网格法）
• FMG（完全多重网格法：嵌套网格与V循环或者W循环结
合）
大规模稀疏矩阵并行计算
4/8/2015
8
代数多重网格法
大规模稀疏矩阵并行计算
4/8/2015
9
代数多重网格法
• 在粗网格上对残差方程进行求解（可用迭代法或直接解
法）。
• 延拓或插值（interpolation）：将细网格节点上的值通过
分片插值延拓到细网格节点上。
• 通过光滑的残差对解进行修正。
• 后光滑（post-smooth），类似于前光滑。
大规模稀疏矩阵并行计算
4/8/2015
10
代数多重网格法方法选择
• 对于非结构化网格形成的矩阵，SGS，SSOR方法不易并
行，即使使用顶点着色技术，因其粗粒度的并行更适合于
传统的多核处理器，并不非常适合GPU这样的细粒度并行
的架构。
• Jacobi方法不具有低通滤波性，因此推荐使用dampJacobi和PCG方法作为迭代子，其中damp-Jacobi方
法的权值一般取为2/3。
• 在最粗网格上的计算推荐使用直接解法。
• 通常对于二阶椭圆边值问题，几何多重网格法具有更好的
计算效率以及收敛速度。
大规模稀疏矩阵并行计算
4/8/2015
11
代数多重网格法方法选择
• 一般遵循两个原则：
o 对于某个顶点，其邻接顶点要么属于粗网格顶点，要么至少连接到一个粗网格顶
点。
o 粗网格顶点集应是任意两个粗网格节点不相邻的极大独立集。
• 有时很难同时满足两个条件，优先满足第一个条件时尽量
满足第二个条件。
大规模稀疏矩阵并行计算
4/8/2015
12
代数多重网格法方法选择
大规模稀疏矩阵并行计算
4/8/2015
13
代数多重网格法的局限性
• 任意几何网格不适用于所有问题。
• 需要高质量的网格划分。
• 不便于编写通用的程序。
• 重点要解决的问题：网格粗化（对应于粗水平方程组）。
• 常用的网格粗化方法复杂：RS，RS2，RS3，Falgout，
HIPS，CLJP。
大规模稀疏矩阵并行计算
4/8/2015
14
大规模稀疏矩阵GPU计算
程序优化设计探索
• 内核执行的优化
o 在大循环中具有大量入口参数的内核，其不变的参数在循环开始前放入常量内存。
避免多余的内存操作
o 合理的网格布局。
o 有时将一个大grid拆分成多个阶段小的grid将有助于提高网格利用率，提高计算效
率，例如对称矩阵的分解以及三角方程组的计算。
• 寄存器优化
o 一个线程中计算输出多个变量，用寄存器内存替换共享内存。
o 在Fermi上，如果程序中存取操作占多数，则对于大于32bit的数据， 以字节流的
形式访问，因为对于例如双精度数据，这时只有一个warp调度器可以工作。
大规模稀疏矩阵并行计算
4/8/2015
15
大规模稀疏矩阵GPU计算
程序优化设计探索
• 合并访问
存取操作以half-warp（计算能力<=1.3的硬件上）或者
warp为单位（Fermi架构上）。
例如，在Fermi架构之前的硬件上，用Vector-CSR方法
计算稀疏矩阵乘法，使用half-warp相对于使用warp有
很高的性能提升。
• 提高线程网格的利用率，减少了线程块的开销。
大规模稀疏矩阵并行计算
4/8/2015
16
大规模稀疏矩阵GPU计算
程序优化设计探索
• 分支优化
消除分支结构的小技巧
例如：
if( a>b ){
a=c;
} else {
a=0;
}
可以替换为：
a=( a>b )*c;
大规模稀疏矩阵并行计算
4/8/2015
17
大规模稀疏矩阵GPU计算
程序优化设计探索
• 指令按照half-warp（计算能力<=1.3）或者warp对齐。
例如：每个线程计算输出7个变量，每个变量的计算差别
很大。这时可以让block的第一个warp的所有线程计算第
一个变量，第二个warp计算第二个变量，……
• 可以利用函数指针（在计算能力<=1.3的硬件上可以使用
对齐到warp边界的控制语句，这时并不会在warp内造成
路径分支（uniform divergence）），通过warp编号来
选择；但是对于相近的计算则不建议使用函数指针反而会
降低效率。
Footer Text
4/8/2015
18
大规模稀疏矩阵GPU计算
程序优化设计探索
• 对于矢量类型数据，使用SOA（Structure of Array）格
式代替
例如，float4可使用 xxxx… yyyy… zzzz… wwww…的存
储结构代替，一般更有效。
• 在Fermi硬件上，读float4类型的数据，虽然显存带宽可
以被充分利用，但是会有部分CUDA Core暂时闲置，并
且必须等待两次的存储请求完成才开始计算，而如果使用
SOA，则在其后的各分量独立的计算中可以更有效隐藏延
迟。
大规模稀疏矩阵并行计算
4/8/2015
19
大规模稀疏矩阵GPU计算
程序优化设计探索
• 如果按照显式的warp模式进行操作，则尽量将每个warp
对应操作的存储器起始地址对齐。如果每个warp的活动
线程数小于75%左右时，则不建议使用。
• 数据结构应该和网格布局相互适应来有效利用存储控制器
的带宽。例如矩阵的转置。
大规模稀疏矩阵并行计算
4/8/2015
20
谢谢！
大规模稀疏矩阵并行计算
4/8/2015
21