這是 標準并行編程 系列的第二篇文章，講述在標準語言中使用并行性來加速計算的優點。

用標準語言并行性開發加速代碼

多個 GPU 標準 C ++并行編程，第 2 部分

將應用程序移植到 GPU 的難度因情況而異。在最佳情況下，您可以通過調用現有的 GPU 優化庫來加速關鍵代碼部分。例如，當模擬軟件的構建塊由 BLAS 線性代數函數組成時，可以使用 cuBLAS 對其進行加速。

但在許多代碼中，你不能四處做一些手工工作。在這些場景中，您可以考慮使用特定于域的語言，例如 CUDA 來針對特定的加速器。或者，您可以使用基于指令的方法，如 OpenMP 或 OpenACC ，以保持原始語言，并使用相同的代碼針對主機和各種類型的設備。

標準并行性

隨著本機形式的并行在C++、FORTRAN和 Python 編程語言的現代版本中的出現，現在可以利用類似的高級方法而不需要語言擴展。

C ++中的標準并行性

我們的重點是 C ++語言，它作為 C ++ 17 標準，在標準庫中提供了許多算法的并行版本。底層編程模型是前面提到的兩種方法的混合體。它的工作方式類似于基于庫的方法，因為 C ++提供了用于排序、搜索和累積和的常見任務的并行算法，并且可以在即將到來的版本中添加對特定領域特定算法的支持。此外，以通用for_each和transform_reduce算法的形式提供了并行手寫循環的形式。

C ++并行算法通過語言的本地語法來代替非標準擴展來表示并行性。通過這種方式，它們保證了所開發軟件的長期兼容性和可移植性。這篇文章還表明，獲得的性能通常與 CUDA C ++等傳統方法獲得的性能相當。

這種方法的新穎之處在于它無縫地集成到現有的代碼庫中。此方法允許您保留軟件體系結構，并有選擇地加快關鍵組件的性能。

將 C ++項目移植到 GPU 可以簡單到通過調用for_each或transform_reduce來替換所有的循環，如果它包含了一個還原。

我們將通過典型的重構步驟來克服當前 C ++編譯器不兼容問題。這篇文章列出了出于性能原因所需的修改，這些修改更具普遍性，原則上獨立于編程形式。這包括以允許合并內存訪問的方式重新構造數據的要求。

對于當前編譯器， C ++并行算法只針對單個 GPU ，而需要明確的 MPI 并行性來針對多個 GPU 。為此，重用現有并行 CPU 代碼的 MPI 后端非常簡單，我們將介紹一些實現最先進性能的指導原則。

我們將討論這些實施規則、指南和最佳實踐，并以 Palabos 為例進行說明， Palabos 是一個基于晶格玻爾茲曼方法（ LBM ）的計算流體力學軟件庫。 PalabOS 在 2021 中被移植到多個 GPU 硬件中，只有幾個月的工作，并且說明了對于原來的代碼很難適應 GPU 的建議重構步驟的需要。由于廣泛使用面向對象的數據結構和編碼機制，原始代碼對 GPU 的適應性較差。

您知道使用 C ++標準并行性允許混合算法，其中一些算法在 GPU 上執行，但有些算法保持在 CPU 上。這完全取決于它們是否適合 GPU 執行，或者僅僅取決于它們的 GPU 端口的進度狀態。這一特性突出顯示了 C ++標準并行性的主要優點之一，與保持原有體系結構和大部分軟件代碼的完整性保持一致。

使用 C ++標準并行編程的 GPU 程序設計

如果這是您第一次聽到 C ++并行算法，您可能想讀 用標準語言并行性開發加速代碼 ，它介紹了 C ++中的標準語言并行性的主題， FORTRAN 和 Python 。

基本概念非常簡單。您可以通過執行策略來獲得許多標準的 C ++算法在主機上或設備上并行運行，作為一個額外的參數提供給算法。在本文中，我們使用par_unseq execution policy，它表示對不同元素的計算是完全獨立的。

以下代碼示例執行并行操作，將std：：vector《double》的所有元素乘以 2 ：

for_each（execution：：par_unseq， begin（v）， end（v）， ［］（double& x） { x *= 2.0;

}）;

該算法由nvc++ compiler和-stdpar option編譯，在 GPU 上執行。根據編譯器、編譯器選項和并行算法的實現，還可以在多核 CPU 或其他類型的加速器上獲得多線程執行。

此示例使用通用的for_each算法，該算法以函數對象的形式將任何元素操作應用于向量 v 。在本例中，它是一個內聯 lambda 表達式。可以使用算法transform_reduce而不是for_each來指定額外的縮減操作。

在for_each算法調用中，調用 lambda 函數時會引用連續的容器元素。但有時，為了訪問外部數據數組或實現非本地模板，還必須知道元素的索引。

這可以通過在 C ++ C ++ 17 中的推力庫（包括NVIDIA HPC SDK ）和std：：ranges：：views：：iota中提供的counting_iterator迭代來完成，或者在 C ++ 20 中更新。在 C ++ 17 中，最簡單的解決方案是從當前元素的地址推導索引。

使用 C ++標準并行性的雅可比示例

為了說明這些概念，下面是一個代碼示例，它使用并行 STL 計算非局部模具操作和誤差估計的縮減操作。它執行雅可比迭代，計算每個矩陣元素的四個最近鄰的平均值：

void jacobi_iteration(vector const& v, vector& tmp) { double const* vptr = v.data(); double *tmp_ptr = tmp.data(); double l2_error = transform_reduce(execution::par_unseq, begin(v), end(v), 0., plus, [=](double& x) { int i = &x - vptr; // Compute index of x from its address. auto [iX, iY] = split(i); double avg = 0.25 * ( vptr[fuse(iX-1, iY)] + vptr[fuse(iX+1, iY)] + vptr[fuse(iX, iY-1)] + vptr[fuse(iX, iY+1)] ); tmp_ptr[i] = avg; return (avg – x) * (avg – x); } );
)

這里，split表示將線性索引i分解為 x 坐標和 y 坐標，fuse則相反。如果域是一個統一的nx-by-ny矩陣， Y 索引在內存中順序運行，則定義如下：

fuse = ［ny］（int iX， int iY） { return iY + ny * iX; }

split = ［ny］（int i） { return make_tuple（i / ny， i % ny）; }

當算法同時執行時，使用臨時向量存儲計算出的平均值可以保證確定性結果。

此代碼的完整和通用版本可從 gitlab.com/unigehpfs/paralg GitLab 存儲庫獲得。該存儲庫還包括一個混合版本（ C ++標準并行和 MPI ），它是圍繞本文提供的建議構建的，在多個 GPU 上高效運行。

您可能已經注意到，沒有明確的語句將數據從主機傳輸到設備，然后再傳輸回來。 C ++標準實際上不提供任何這樣的語句，并且任何并行算法在一個設備上的實現必須依賴于自動存儲器傳輸。通過NVIDIA HPC SDK ，這是通過 CUDA 統一內存，從 CPU 和 GPU 訪問的單個存儲器地址空間來實現的。如果代碼在 CPU 上訪問此地址空間中的數據，然后在 GPU 上訪問，則內存頁會自動遷移到訪問處理器。

對于使用 GPU 加速 CPU 應用程序， CUDA 統一內存特別有用，因為它使您能夠專注于以增量方式逐個函數移植應用程序的算法，而無需擔心內存管理。

另一方面，隱藏數據傳輸的性能開銷很容易抵消 GPU 的性能優勢。通常，在 GPU 上生成的數據應盡可能保存在 GPU 內存中，通過并行算法調用表示其所有操作。這包括數據后處理，如數據統計計算和可視化。如本文 Part 2 所示，它還包括 MPI 通信的數據打包和解包。

按照這篇文章的建議，將代碼移植到 GPU 變得非常簡單，只需通過調用并行算法來替換所有時間關鍵的循環和相關的數據訪問。不過，最好記住， GPU 通常比 CPU 擁有更多的內核，并且應該暴露在更高級別的并行性中。例如，在下一節介紹的流體動力學問題中，流體域被均勻的、類似矩陣的網格部分覆蓋。

在原始的 CPU 代碼中，每個 CPU 核心按順序處理一個或多個網格部分，如圖 1 頂部所示。至于 GPU ，網格部分元素上的環路應該并行，以完全占據 GPU 核心。

示例： Lattice Boltzmann 軟件和 Palabos

LBM 采用顯式時間步格式求解流體流動方程，涵蓋了廣泛的應用。這包括經過復雜幾何形狀的流動，如多孔介質、多相流、可壓縮超音速流等。

LBM 通常比其他解算器在數值網格的每個節點上分配更多變量。當經典的不可壓縮 Navier-Stokes 解算器僅用三個變量表示速度分量，外加一個臨時壓力項時， LBM 方法通常需要 19 個變量，稱為populations。因此， LBM 的內存占用空間要高出 5-6 倍。如果 Navier Stokes 解算器使用臨時變量，或者如果系統中添加了進一步的物理量（如密度和溫度），實際里程可能會有所不同。

因此，豐富的內存訪問和較低的運算強度是 LBM 的特點。在集群級 GPU ，像 NVIDIA V100 和 NVIDIA A100 一樣，性能完全受限于內存訪問，甚至對于計算密集和復雜的 LBM 方案也是如此。

以 NVIDIA A100 40 GB GPU為例，它具有1555 Gb/s的內存帶寬。在每一個明確的時間步長，每個節點訪問19個變量或［EZX27 ］，每個都占用八個字節，每個都是雙精度的。它們被計數兩次：一次用于從 GPU 內存到 GPU 內核的數據傳輸，另一次用于在計算操作后寫回 GPU 內存。

假設一個完美的內存子系統和最大的數據重用， LBM 的峰值吞吐量性能為每秒處理 1555 /（ 19 * 8 * 2 ）= 51.1 億個網格節點。在 LBM 術語中，通常使用每秒千兆晶格節點更新（ GLUPS ），例如 5.11 GLUPS 。

然而，在現實生活中的應用程序中，每個節點都會額外讀取一些信息來管理域異常情況。在 Palabos 中，這是節點標記的 32 位整數和額外數據數組的 64 位索引，有效地將峰值性能降低到 4.92 GLUPS 。

該模型提供了一種簡單的方法來估計 LBM 代碼可以達到的最佳峰值性能，因為緩存中不適合足夠大的網格。我們在整個帖子中使用這個模型來證明用 C ++并行算法獲得的性能是一樣好的。在幾個百分點的差距之外，無論是 CUDA 、 OpenMP ，還是任何其他 GPU 形式主義，都不能做得更好。

LBM 巧妙地區分了由局部collision step表示的計算和封裝在streaming step中的內存傳輸操作。以下代碼示例顯示了具有矩陣式拓撲結構的結構化網格的典型時間迭代：

for (int i = 0; i < N; ++i) { // Fetch local populations "f" from memory. double f_local[19]; for (int k = 0; k < 19; ++k) { f_local[k] = f[i][k]; } collide(f_local); // Execute collision step. // Write data back to neighboring nodes in memory (streaming). auto [iX, iY, iZ] = split(i); for (int k = 0; k < 19; ++k) { int nb = fuse(iX+c[k][0], iY+c[k][1], iZ+c[k][2]); ftmp[nb][k] = f_local[k]; }
}

與前一節中的 Jacobi 迭代一樣，該函數將計算出的數據寫入temporary array ftmp，以避免多線程執行期間出現爭用情況，這使其成為演示本文概念的理想候選。有 替代就地算法 可以避免內存復制。然而，它們更復雜，因此不太適合用于說明目的。

自然過程的模擬和建模 課程介紹了LBM。有關如何使用C++并行算法開發 GPU 的LBM代碼的更多信息，請參見多核格子玻爾茲曼模擬的跨平臺編程模型。

在本文中，我們使用 開源 LBM 庫 Palabos 來展示如何將現有的 C ++庫用并行算法移植到多 GPU 。乍一看， Palabos 似乎不適合 GPU 端口，因為它強烈依賴面向對象的機制。然而，在 Palabos 的案例中，我們將介紹幾種變通方法，這些方法只需對代碼體系結構進行表面的更改即可實現最先進的性能。

從面向對象設計轉向面向數據設計

為了服務于大型社區，Palabos強調多態性和其他面向對象技術。包含數據（種群）和方法（局部碰撞模型）的對象代表每個網格節點。這為開發新模型提供了一個方便的API，并提供了一個靈活的機制來調整模型從一個單元到另一個單元的物理行為或數值方面。

然而，由于數據布局效率低下、執行路徑復雜，以及對虛擬函數調用的依賴，這種面向對象的方法不太適合在 GPU 上執行。以下幾節將教您如何通過采用開發模型，以 GPU 友好的方式重構代碼，我們在總括術語下稱之為data-oriented programming。

擺脫基于類的多態性

圖 2 的左半部分展示了 Palabos 中網格節點上沖突模型的典型代碼執行鏈。算法的不同組件被疊加起來，并通過虛擬函數調用，包括底層的數值 LBM 算法（ RR ）、附加物理（“ Smagorinsky ”）和附加數值方面（左邊界）。

這個面向對象設計的教科書案例變成了代碼 GPU 端口的責任。這個問題的出現是因為當前版本的 HPCSDK 不支持 C ++并行算法中的 GPU 的虛函數調用。一般來說，出于性能原因，在 GPU 上應該避免這種類型的設計，因為它限制了執行路徑的可預測性。

簡單的解決方法是將執行鏈收集到單個函數中，這些函數按順序顯式調用各個組件，并用唯一的標記標識它們。

在 Palabos 中，這個獨特的標記是在序列化機制的幫助下生成的，該機制最初是為了支持動態自適應仿真的檢查點和網絡通信而開發的。這表明，如果重構軟件項目的體系結構足夠靈活，那么 GPU 端口的大部分重構工作都是自動完成的。

現在，您可以為每個網格節點提供一個標記，標識完整沖突步驟的代碼，并用一個大的 switch 語句表示沖突步驟：

switch（tag） { case rr_les： fun_rr_les（f_local）; break; case rr_les_BCleft： fun_rr_les_BCleft（f_local）; break; …

}

隨著 switch 語句變大，由于生成的內核在 GPU 內存中占用的空間，它可能會遇到性能問題。

另一個問題是軟件項目的可維護性。目前，如果不修改這個 switch 語句（它是庫核心的一部分），就不可能提供新的碰撞模型。這兩個問題的解決方案在于，在編譯時使用 C ++模板機制在最終用戶應用程序中生成具有選定數量的實例的切換語句。 Palabos-GPU 資源頁詳細介紹了這種技術。

重新安排內存以鼓勵合并的內存訪問

面向對象的設計還會導致內存布局無法在 GPU 的多核架構上高效處理。當每個節點對 LBM 方案的 19 個局部種群進行分組時，數據以結構陣列（ AoS ）的形式結束（圖 3 ）。為了簡單起見，每個節點只顯示四個總體。

AoS 數據布局導致性能不佳，因為它阻止了流式處理步驟中的合并內存訪問，而流式處理步驟由于非本地模板，在內存訪問方面是算法最關鍵的部分。

數據應該以陣列結構（ SoA ）的方式對齊，在流式處理步驟中進行通信的給定類型的所有群體都在連續的內存地址上對齊。在這種重新安排之后，即使是一個相對簡單的 LBM 算法也能獲得接近典型 GPU 內存帶寬的 80% 。

面向數據的設計意味著您將中心重要性賦予數據的結構和布局，并圍繞此結構構建數據處理算法。面向對象的方法通常采用反向路徑。

GPU 端口的第一步應該是了解應用程序的理想數據布局。就 LBM 而言， GPU 上 SoA 布局的優越性是眾所周知的事實。內存布局和內存遍歷算法的細節在之前發布的案例研究 開源 STLBM 代碼 中進行了測試。

結論

在這篇文章中，我們討論了使用 C ++標準并行編程編寫 GPU 應用程序的基本技術。我們還提供了晶格玻爾茲曼方法和 Palabos 應用程序的背景信息，我們在案例研究中使用了這些信息。最后，我們討論了兩種方法，可以重構源代碼，使其更適合在 GPU 上運行。

在下一篇文章中，我們繼續使用這個應用程序，并討論如何在NVIDIA上運行時在C++應用程序中獲得高性能。我們還演示了如何通過MPI擴展應用程序以使用多個 GPU 。

關于作者

Jonas Latt 是瑞士日內瓦大學計算機科學系的副教授。他從事高性能計算和計算流體力學的研究，并在包括地球物理、生物醫學和航空航天領域在內的跨學科領域進行應用。他是 lattice Boltzmann 復雜流動模擬開源軟件 Palabos 的最初開發者和當前共同維護者。他以前在日內瓦大學獲得物理學和計算機科學博士學位，并通過塔夫斯大學（波士頓，美國）和綜合理工學校 F.E.EdRaelde 洛桑 EPFL （瑞士）的研究，并作為 CFD 公司 FuluKIT 的聯合創始人，對流體力學感興趣。

Christophe Guy Coreixas 是一名航空工程師， 2014 年畢業于 ISAE-SUPAERO （法國圖盧茲）。 2018 年，他在 CERFACS 從事面向行業應用的可壓縮晶格玻爾茲曼方法研究時獲得了博士學位（流體動力學）。作為日內瓦大學計算機科學系的博士后，克里斯多夫現在開發格子玻爾茲曼模型來模擬航空、多物理和生物醫學流程。

Gonzalo Brito 是 NVIDIA 計算性能與 HPC 團隊的高級開發技術工程師，工作于硬件和軟件的交叉點。他熱衷于讓加速計算變得更容易實現。在加入NVIDIA 之前，岡薩洛在 RWTH 亞琛大學空氣動力學研究所開發了多物理方法，用于顆粒流。

Jeff Larkin 是 NVIDIA HPC 軟件團隊的首席 HPC 應用程序架構師。他熱衷于高性能計算并行編程模型的發展和采用。他曾是 NVIDIA 開發人員技術小組的成員，專門從事高性能計算應用程序的性能分析和優化。 Jeff 還是 OpenACC 技術委員會主席，曾在 OpenACC 和 OpenMP 標準機構工作。在加入NVIDIA 之前，杰夫在位于橡樹嶺國家實驗室的克雷超級計算卓越中心工作。

審核編輯：郭婷