大規模稀疏矩陣并行計算

李修宇QQ：2955533816/29/20231主流求解方法直接法GAUSS消去法波前法多波前法迭代法經典迭代法Jacobi、SOR、SSOR投影方法CG、GMRES預處理技術不完全分解預處理條件代數多重網格技術6/29/2023大規模稀疏矩陣并行計算2矩陣性質對求解的影響性質影響6/29/2023大規模稀疏矩陣并行計算3非零元的分布帶狀分布按塊分布……正定性對稱性矩陣的存儲方式求解方法的選擇求解速度……直接法矩陣圖重排：一般分為兩大類，帶寬縮減算法（也常稱為外形縮減）和區域分解算法，應用較多的帶寬縮減算法CM，RCM，GPS，Rosen算法。一般建議多重方法結合使用：全局方法的全局平衡性、局部方法的局部最優特性。符號分解：確定非零元結構以及相應的消元索引，以便在實際數值分解前確定所需存儲資源大小，避免數值分解中動態分配存儲空間和復雜的索引策略。構建消去樹(eliminationtree)：確定分解節點之間的分解依賴，即確定分解的順序并構成并行分解的層次結構。6/29/2023大規模稀疏矩陣并行計算4

直接法數值分解：利用符號分解得到的非零元結構和索引沿消去樹路徑進行分解?；卮蠼猓喊ㄇ跋颍╢orward）和后向（backward）回代，可先構建消去依賴樹或頂點著色技術實現并行回代求解。在有限元領域應用最廣的直接求解方法常使用帶寬縮減或多區域分解的多波前法（multifrontal）。6/29/2023大規模稀疏矩陣并行計算5對稱正定矩陣的求解

6/29/2023大規模稀疏矩陣并行計算6對稱矩陣的不完全分解

6/29/2023大規模稀疏矩陣并行計算7代數多重網格法V-CycleAMG（V循環多重網格法）W-CycleAMG（W循環多重網格法）FMG（完全多重網格法：嵌套網格與V循環或者W循環結合）6/29/2023大規模稀疏矩陣并行計算8代數多重網格法

6/29/2023大規模稀疏矩陣并行計算9代數多重網格法在粗網格上對殘差方程進行求解（可用迭代法或直接解法）。延拓或插值（interpolation）：將細網格節點上的值通過分片插值延拓到細網格節點上。通過光滑的殘差對解進行修正。后光滑（post-smooth），類似于前光滑。6/29/2023大規模稀疏矩陣并行計算10代數多重網格法方法選擇對于非結構化網格形成的矩陣，SGS，SSOR方法不易并行，即使使用頂點著色技術，因其粗粒度的并行更適合于傳統的多核處理器，并不非常適合GPU這樣的細粒度并行的架構。Jacobi方法不具有低通濾波性，因此推薦使用damp-Jacobi和PCG方法作為迭代子，其中damp-Jacobi方法的權值一般取為2/3。在最粗網格上的計算推薦使用直接解法。通常對于二階橢圓邊值問題，幾何多重網格法具有更好的計算效率以及收斂速度。6/29/2023大規模稀疏矩陣并行計算11代數多重網格法方法選擇一般遵循兩個原則：對于某個頂點，其鄰接頂點要么屬于粗網格頂點，要么至少連接到一個粗網格頂點。粗網格頂點集應是任意兩個粗網格節點不相鄰的極大獨立集。有時很難同時滿足兩個條件，優先滿足第一個條件時盡量滿足第二個條件。6/29/2023大規模稀疏矩陣并行計算12代數多重網格法方法選擇

6/29/2023大規模稀疏矩陣并行計算13代數多重網格法的局限性任意幾何網格不適用于所有問題。需要高質量的網格劃分。不便于編寫通用的程序。重點要解決的問題：網格粗化（對應于粗水平方程組）。常用的網格粗化方法復雜：RS，RS2，RS3，Falgout，HIPS，CLJP。6/29/2023大規模稀疏矩陣并行計算14大規模稀疏矩陣GPU計算程序優化設計探索內核執行的優化在大循環中具有大量入口參數的內核，其不變的參數在循環開始前放入常量內存。避免多余的內存操作合理的網格布局。有時將一個大grid拆分成多個階段小的grid將有助于提高網格利用率，提高計算效率，例如對稱矩陣的分解以及三角方程組的計算。寄存器優化一個線程中計算輸出多個變量，用寄存器內存替換共享內存。在Fermi上，如果程序中存取操作占多數，則對于大于32bit的數據，以字節流的形式訪問，因為對于例如雙精度數據，這時只有一個warp調度器可以工作。6/29/2023大規模稀疏矩陣并行計算15大規模稀疏矩陣GPU計算程序優化設計探索合并訪問存取操作以half-warp（計算能力<=1.3的硬件上）或者warp為單位（Fermi架構上）。例如，在Fermi架構之前的硬件上，用Vector-CSR方法計算稀疏矩陣乘法，使用half-warp相對于使用warp有很高的性能提升。提高線程網格的利用率，減少了線程塊的開銷。6/29/2023大規模稀疏矩陣并行計算16大規模稀疏矩陣GPU計算程序優化設計探索分支優化消除分支結構的小技巧

例如： if(a>b){ a=c; }else{ a=0; }可以替換為： a=(a>b)*c;6/29/2023大規模稀疏矩陣并行計算17大規模稀疏矩陣GPU計算程序優化設計探索指令按照half-warp（計算能力<=1.3）或者warp對齊。例如：每個線程計算輸出7個變量，每個變量的計算差別很大。這時可以讓block的第一個warp的所有線程計算第一個變量，第二個warp計算第二個變量，……可以利用函數指針（在計算能力<=1.3的硬件上可以使用對齊到warp邊界的控制語句，這時并不會在warp內造成路徑分支（uniformdivergence）），通過warp編號來選擇；但是對于相近的計算則不建議使用函數指針反而會降低效率。6/29/2023FooterText18大規模稀疏矩陣GPU計算程序優化設計探索對于矢量類型數據，使用SOA（StructureofArray）格式代替例如，float4可使用xxxx…yyyy…zzzz…wwww…的存儲結構代替，一般更有效。在Fermi硬件上，讀float4類型的數據，雖然顯存帶寬可以被充分利用，但是會有部分CUDACore暫時閑置，并且必須等待兩次的存儲請求完成才開始計算，而如果使用SOA，則在其后的各分量獨立的計算中可以更有效隱藏延遲。6/29/2023大規模稀疏矩陣并行計算19大規模稀疏矩陣GPU計算程序優化設計探索如果按照顯式的warp模式進行操