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