これは「標準並列プログラミング」シリーズの 2 回目の投稿です。このシリーズでは、標準言語による並列化をアクセラレーテッド コンピューティングで使用することの利点を取り上げています。
アプリケーションを GPU に移植することの難しさはケースバイケースです。最高のシナリオは、GPU で最適化された既存のライブラリを呼び出すことで、重要なコード セクションを高速化することです。たとえば、シミュレーション ソフトウェアの構成要素が BLAS の線形代数関数で構成されている場合、cuBLAS を利用して高速化することが可能です。
しかしながら、多くのコードでは、相当な量の手作業を回避できません。そのような場合は、特定のアクセラレータをターゲットにした CUDA のようなドメイン固有言語の使用を検討するでしょう。あるいは、OpenMP や OpenACC などのディレクティブベースの手法を利用することで元の言語のまま、同じコードで、ホストとさまざまな種類のデバイスをターゲットにできます。
標準並列処理
C++、Fortran、Python などのプログラミング言語の最新版でネイティブ形式の並列処理が登場したことで、言語拡張に頼ることなく、同様のハイレベルな手法を活用できるようになりました。
C++ における標準並列処理
重点的に取り組んでいるものが C++ です。この言語では、C++17 以降、標準ライブラリでさまざまなアルゴリズムの並列バージョンを提供しています。基礎となるプログラミング モデルは、先に述べた 2 つの手法を組み合わせたものです。C++ は並べ替え、検索、累積和などの一般的なタスクに並列アルゴリズムを提供し、また今後のバージョンでドメイン固有のアルゴリズムを追加する可能性があります。そのような理由から、部分的にはライブラリベースのアプローチのように機能します。さらに、並列の手書きループは、一般的な for_each アルゴリズムと transform_reduce アルゴリズムの形式となります。
C++ 並列アルゴリズムは、非標準の拡張機能の代わりに、言語のネイティブ構文を介して並列処理を表現します。そのため、開発したソフトウェアは、長期間にわたって互換性と移植性が保証されます。この投稿では、得られたパフォーマンスが一般的に CUDA C++ などの従来の手法で得られたパフォーマンスに匹敵することも示します。
この手法の斬新な点は、既存のコード ベースにシームレスに統合されることです。この方法により、ソフトウェア アーキテクチャを維持したまま、重要なコンポーネントのパフォーマンスを選択的に高速化できます。
C++ プロジェクトを GPU に移植することは、全ての for ループを for_each、リダクションが含まれる場合は transform_reduce の呼び出しに置換することと同じくらい簡単です。
現行の C++ コンパイラの非互換性問題を克服するために必要な典型的リファクタリング手順を説明します。この投稿には、パフォーマンス上の理由で必要な一連の変更が含まれています。それらの変更は汎用性の高いものであり、原則として、プログラミング形式には依存しません。これには、コアレス メモリ アクセスを可能にする方法でデータを再構築するための要件が含まれています。
現行のコンパイラでは、C++ 並列アルゴリズムは単一 GPU のみをターゲットにします。複数の GPU をターゲットにするには、明示的な MPI 並列が必要になります。その目的のためには、既存の並列 CPU コードの MPI バックエンドを再利用するのが簡単です。最先端のパフォーマンスを得るためのガイドラインをいくつか提示します。
実装のルール、ガイドライン、ベスト プラクティスからなるこのリストについて、Palabos の例で説明します。これは、LBM (格子ボルツマン法) に基づく数値流体力学のソフトウェア ライブラリです。Palabos は 2021 年にわずか数か月の作業でマルチ GPU ハードウェアに移植されました。元のコードが GPU にうまく適していないため、推奨リファクタリング手順の必要性が示されています。元のコードが GPU にうまく適していないのは、オブジェクト指向のデータ構造とコーディング メカニズムが広く使用されていることにあります。
C++ 標準並列処理を使用すると、一部のアルゴリズムを GPU で実行し、他を CPU に留めるハイブリッド手法が可能になります。これはすべて、GPU 実行に適しているかに依存します。あるいは単純に、GPU 移植の進行状況に依存します。この特徴は、C++ 標準並列処理の大きな利点のひとつとして際立っています。元のソフトウェア プロジェクトのグローバル アーキテクチャとほとんどのコード ベースをそのまま維持できるようなものです。
C++ 標準並列処理を利用した GPU プログラミング
C++ 並列アルゴリズムについて初めて耳にする場合、「標準言語による並列処理を利用したコードの高速化」をお読みください。C++、Fortran、Python などの標準言語による並列処理を取り上げています。
基本的なコンセプトは驚くほどシンプルです。多くの標準 C++ アルゴリズムは、実行ポリシーを利用することで、ホストまたはデバイスで並列実行できます。実行ポリシーは、アルゴリズムの追加の引数として提供されます。この投稿では、par_unseq 実行ポリシーを使用しています。これは、各要素での計算が完全に独立であることが表現されています。
次のコード例では、std::vector<double> の全要素を二倍にする並列演算が実行されます。
for_each(execution::par_unseq, begin(v), end(v), [](double& x) {
x *= 2.0;
});nvc++ コンパイラと -stdpar オプションでコンパイルされたこのアルゴリズムは GPU で実行されます。コンパイラ、コンパイラ オプション、並列アルゴリズムの実装に応じて、マルチコア CPU またはその他の種類のアクセラレータでマルチスレッド実行することもできます。
この例では、関数オブジェクトの形式でベクトル v に要素ごとの演算を適用する一般的な for_each アルゴリズムが使用されます。この場合、インライン ラムダ式になります。追加のリダクション演算は、アルゴリズム transform_reduce を for_each の代わりに使用することで指定できます。
for_each アルゴリズムの呼び出しでは、連続するコンテナー要素を参照してラムダ関数が呼び出されます。ただし、外部データ配列にアクセスするためや非局所なステンシルを実装するために、要素のインデックスを知らなければならないこともあります。
これは、C++17 の Thrust ライブラリ (NVIDIA HPC SDK に付属) と C++20 以降の std::ranges::views::iota で提供される counting_iterator で反復処理することで行えます。C++17 の場合、最も単純な解決策は、現在の要素のアドレスからインデックス推測することです。
C++ 標準並列処理を利用した Jacobi の例
こうした概念を説明するためのコード例があります。並列 STL を使用し、非局所的な ステンシル演算とエラー見積もりのリダクション演算を実行します。Jacobi 反復法を実行し、すべての行列要素の最近接 4 個の平均値が評価されます。
void jacobi_iteration(vector<double> const& v, vector<double>& 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>, [=](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 を使用すると、これは、CPU と GPU の両方からアクセスできる単一のメモリ アドレス空間である CUDA ユニファイド メモリで達成されます。コードが CPU 上でこのアドレス空間のデータにアクセスし、次に GPU 上でアクセスする場合、アクセスするプロセッサにメモリ ページが自動的に転送されます。
特に GPU を利用して CPU アプリケーションを高速化する場合に、CUDA ユニファイド メモリは便利です。アプリケーションのアルゴリズムをインクリメンタルに、関数ごとに、メモリ管理を心配せずに移植できるためです。
マイナス面は、隠されたデータ転送が余計な負荷を生み、GPU のメリットを相殺してしまうという点です。原則として、GPU で生成されるデータは、可能な限り、GPU メモリに保持するべきです。その方法として、並列アルゴリズム呼び出しでその操作をすべて表現します。これには、データ統計の計算や可視化など、データの後処理が含まれます。この投稿のパート 2 にあるように、MPI 通信のデータのパッキングとアンパッキングも含まれます。
この投稿の推奨事項に従えば、時間がかかるループや関連データ アクセスを並列アルゴリズムの呼び出しに置き換えるだけで、簡単にコードを GPU に移植することができます。ただし、GPU は通常、CPU よりもはるかに多くのコアを処理します。GPU は、より上のレベルの並列処理に利用すべきです。次のセクションで提示される流体力学の問題では、たとえば、均一で行列のようなグリッド部分によって流体ドメインが覆われています (図 1)。
元の CPU コードでは、各 CPU コアが 1 つまたは複数のグリッド部分を順番に処理します。これは図 1 の上部に示されています。GPU については、グリッド部分の要素のループを並列化し、GPU コアを完全に占有する必要があります。
例: 格子ボルツマン ソフトウェアと Palabos
LBM は、明示的な時間ステップ方式で流体の流れの方程式を解き、幅広い応用範囲をカバーします。これには、多孔質媒体、混相流、圧縮性超音速流など、複雑な形状を通過する流れが含まれます。
LBM では多くの場合、数値メッシュのすべてのノードに他のソルバーよりも多くの変数を割り当てます。従来の非圧縮性の Navier-Stokes ソルバーが、速度成分のわずか 3 つの変数と一時的な圧力項の 1 つの変数で状態を表すのに対し、LBM 手法では通常、ポピュレーションと呼ばれる 19 の変数が必要になります。そのため、LBM のメモリ フットプリントは 5 倍から 6 倍になります。Navier-Stokes ソルバーで一時的変数を使用する場合や、密度や温度など、さらなる物理量が追加される場合、実際の量は変わることがあります。
その結果、豊富なメモリ アクセスと低い演算強度が LBM の特徴となります。NVIDIA V100 や NVIDIA A100 など、クラスターレベルの GPU では、計算負荷が高く、洗練された LBM スキームであっても、パフォーマンスはメモリ アクセスによって完全に制限されます。
たとえば、NVIDIA A100 40 GB GPU のメモリ帯域幅は 1555 GB/s です。明示的なタイム ステップごとに、各ノードは 19 個の変数ポピュレーションにアクセスし、それぞれが倍精度で 8 バイトを占有します。これらは 2 回カウントされます。GPU メモリから GPU コアにデータを転送するために 1 回、演算後に再度、GPU メモリに書き込むために 1 回です。
完璧なメモリ サブシステムと最大のデータ再利用を前提とすると、LBM のピーク スループット パフォーマンスでは 1555/(19 * 8 * 2) = 毎秒 51 億 1000 万のグリッド ノードを処理します。LBM では一般的な習慣として、5.11 GLUPS のように、Giga Lattice-node Update Per Second (GLUPS) を使用します。
しかしながら、実際のアプリケーションでは、各ノードによって付加的にいくつかの情報が読み込まれ、ドメインの不規則性に対処します。Palabos では、これはノード タグの 32 ビット整数であり、余分なデータ配列の 64 ビット インデックスであり、実質的にピーク パフォーマンスを 4.92 GLUPS まで減らします。
十分に大きなグリッドがキャッシュに収まらないため、このモデルは LBM コードで達成しうるベスト ピーク パフォーマンスを見積もる簡単な方法を与えます。このモデルは本投稿全体で使用し、C++ 並列アルゴリズムで得られるパフォーマンスが可能な限り良いパフォーマンスであるということを実証します。CUDA も、OpenMP も、その他の GPU フォーマリズムも、数パーセントの差以上のものを与えることはありません。
LBM は、ローカルな衝突ステップで表される計算と、並進ステップでカプセル化されるメモリ転送操作を正確に区別します。次のコード例は、行列状のトポロジを持つ構造化メッシュの場合の典型的な時間反復を示しています。
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 反復法と同様に、この関数は計算後のデータを一時的な配列 ftmp に書き込み、マルチスレッド実行中の競合状態を回避します。この投稿の概念を実証するものとして最適です。メモリの重複を回避する代替のインプレース アルゴリズムがあります。しかしながら、ずっと複雑になるため、説明の目的には適していません。
LBM 入門は「Simulation and modeling of natural processes」コースにあります。C++ 並列アルゴリズムを使用して GPU 用 LBM コードを開発する方法については、「Cross-platform programming model for many-core lattice Boltzmann simulations」を参照してください。
この投稿では、オープンソース LBM ライブラリ Palabos を使用し、既存の C++ ライブラリを並列アルゴリズムでマルチ GPU に移植する方法を紹介します。一見すると、Palabos はオブジェクト指向のメカニズムに強く依存しているため、GPU 移植には適さないように思われます。そこで、回避策をいくつか紹介します。Palabos の場合、コード アーキテクチャを表面的に変更するだけで最先端のパフォーマンスが得られます。
オブジェクト指向からデータ指向の設計に切り替える
大規模なコミュニティに対応するため、Palabos では、ポリモーフィズムとその他のオブジェクト指向手法が重視されています。データ (ポピュレーション) とメソッド (ローカルな衝突モデル) の両方を保持するオブジェクトは、あらゆるグリッド ノードを表します。これにより、新しいモデルを開発するための便利な API と、モデルの物理的動作や数値的側面をセル間で調整するための柔軟なメカニズムを提供します。
しかしながら、このオブジェクト指向の手法は、非効率なデータ レイアウト、複雑な実行パス、仮想関数呼び出しへの依存に起因し、GPU での実行にはあまり適していません。次のセクションでは、GPU フレンドリーな方法でコードをリファクタリングする方法を説明します。具体的には、データ指向プログラミングと呼んでいる開発モデルを採用します。
クラスベースのポリモーフィズムを取り除く
図 2 の左側は、Palabos のグリッド ノードの衝突モデルに対する典型的なコード実行チェーンを示しています。アルゴリズムのさまざまなコンポーネントが積み重ねられ、基礎となる数値 LBM アルゴリズム (RR)、追加の物理的性質 (「Smagorinsky」)、追加の数値的側面 (左境界) など、仮想関数呼び出しによって呼び出されます。
この標準的なオブジェクト指向設計は、コードの GPU 移植にとってやっかいな問題です。この問題は、現行版の HPC SDK で 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++ テンプレート メカニズムを使用して、コンパイル時にエンドユーザー アプリケーションで選ばれたケース数を持つ switch ステートメントを生成することです。この手法の詳細は「Palabos-GPU」リソース ページにあります。
コアレス メモリ アクセスを増やすためのメモリ再配置
オブジェクト指向設計はまた、GPU の多コア アーキテクチャでは効率的に処理できないメモリ レイアウトにつながります。すべてのノードで LBM スキームの 19 個のローカル ポピュレーションをグループ化すると、データは Array-of-Structure (AoS) 配置になります (図 3)。わかりやすくするために、ノードあたり 4 つのポピュレーションのみを示しています。
AoS (Array-of-Structure) データ レイアウトの場合、並進ステップ中にコアレス メモリ アクセスが防止されるため、パフォーマンスが低下します。並進ステップは非局所的なステンシルのため、メモリ アクセスに関してアルゴリズムの最重要部分になります。
データは、並進ステップ中に通信する特定の種類の全ポピュレーションが、連続するメモリ アドレスで整列される SoA (Structure-of-Arrays) 方式で整列されるべきです。この再配列の後、比較的単純な LBM アルゴリズムであっても、典型的な GPU のメモリ帯域幅の 80% に近くなりました。
データ指向設計は、データの構造とレイアウトを重視し、データ構造を中心にデータ処理アルゴリズムを構築することを意味します。オブジェクト指向の手法では通常、方針が逆になります。
GPU 移植の最初の手順は、アプリケーションにとって理想的なデータ レイアウトを知ることになります。LBM の場合、GPU では SoA レイアウトのほうが優れたレイアウトになります。オープンソース STLBM コードとして公開されている先行事例では、メモリ レイアウトとメモリ トラバーサル アルゴリズムの仕様がテストされました。
結論
この投稿では、C++ 標準並列プログラミングを使用した GPU アプリケーション開発の基本的な手法について説明しました。また、格子ボルツマン法と、研究事例で使用した Palabos アプリケーションの背景を説明しました。最後に、GPU 実行に適するように、ソース コードをリファクタリングする 2 通りの方法について説明しました。
次の投稿では、引き続きこのアプリケーションを取り上げ、NVIDIA GPU で実行するとき、C++ アプリケーションで高いパフォーマンスを得る方法について説明します。また、MPI で複数の GPU を使用するようにアプリケーションを拡張する方法について紹介します。
