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