はじめに

記事「MPIによる分散通信を試す」ではMPIの簡単な使い方について説明しました。
今回はOpenACCと組み合わせて、マルチノードかつマルチGPUによる並列計算の実装を紹介します。

MPIは1ノード内で複数プロセスを起動でき、かつ1プロセスだけの実行も可能なので、OpenACC/MPIは以下を1つの実装で実現可能です。
（無条件ではなく、そのようにコードが実装されている必要があります）

• シングルノード、シングルGPU
• シングルノード、マルチGPU
• マルチノード、マルチGPU

MPI上でのGPUデバイス割り当て

OpenACCではacc_set_device_num(device_number, acc_device_nvidia)関数で利用するGPUを切り替え、カーネル実行やメモリコピーを制御していました。
OpenACC/MPIでもacc_set_device_num、GPU数を返すacc_get_num_devices(acc_device_nvidia)関数と、MPIプロセス数とランク番号を組み合わせることでプロセスごとに利用するGPUを切り替えることが可能です。

OpenACCのAPIは1ノード上のGPUに対する操作を提供しているため、acc_get_num_devicesは各MPIプロセスが実行されているノードに存在するGPUの数を返すことに注意してください。
一方、MPIプロセスのランク番号はmpirunで実行したすべてのプロセスへの連番なので、グローバルなランク番号とノードローカルなGPU番号を対応付けする必要があります。
各ノードに接続されているGPUの数が同一であることを仮定すると、以下のような実装になります。

#include <mpi.h>
#include <openacc.h>
#include <assert.h>
#include <stdio.h>
int main(int argc, char *argv[]) {
    MPI_Init(&argc, &argv);
    int rank, nprocs;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
    const int ngpu_per_node = acc_get_num_devices(acc_device_nvidia);
    assert(ngpu_per_node > 0);
    // gpu_id = [0, 1, 2, ..., ngpu_per_node - 1]
    const int gpu_id = rank % ngpu_per_node;
    acc_set_device_num(gpu_id, acc_device_nvidia);
    printf("Rank %d/%d = GPU %d/%d\n", rank, nprocs, gpu_id, ngpu_per_node);
    MPI_Finalize();
}

NVIDIA HPC SDKのOpen MPI/HPC-Xのコンパイラを利用して、OpenACCが有効な状態でMPIプログラムをコンパイルします。
オプションは通常のOpenACCプログラムと同様です。

$ mpicxx -acc -o ./xtest acc_mpi_multi_gpu_0.cc

mpirunはデフォルトでは実行したノードで全プロセスを起動しますが、ホスト名と実行するプロセス数を指定したホストファイルを渡すことで、別ノードで実行することが可能です。
ホスト名の設定方法はシステム環境により異なりますが、仮にhost-XX.example.comとします。XXにはノード番号が入ります。
ホストファイルの書き方もMPIコンパイラやmpirunと同様、処理系に応じて異なります。Open MPI/HPC-Xの場合はホスト名 + slots=Nを1行ずつ記載します。

例えば2ノードで各ノード4プロセスを実行する場合、以下のホストファイルを作成します。
このプロセス数はOpenACC/MPIの場合、ノードにあるGPU数に揃えます。つまり2ノード4 GPUの合計8 GPUで実行することにします。

$ cat hostfile.txt
host-00.example.com slots=4
host-01.example.com slots=4

このホストファイルを使って8プロセスで実行すると、先のOpenACC/MPIプログラムの実行ログは以下のようになります。

$ mpirun -n 8 --hostfile hostfile.txt --bind-to none ./xtest | sort
Rank 0/8 = GPU 0/4
Rank 1/8 = GPU 1/4
Rank 2/8 = GPU 2/4
Rank 3/8 = GPU 3/4
Rank 4/8 = GPU 0/4
Rank 5/8 = GPU 1/4
Rank 6/8 = GPU 2/4
Rank 7/8 = GPU 3/4

少し古いGPU対応のMPI処理系では、MPI_Initを呼ぶ前にプロセスが利用するGPU割り当てが必要となる場合があります。
Open MPI/HPC-Xでは、mpirunで実行したプログラムにノードローカルなランク番号とランクサイズが記載された2つの環境変数が定義されています。

• OMPI_COMM_WORLD_LOCAL_RANK
• OMPI_COMM_WORLD_LOCAL_SIZE

先ほどのコードと同等処理をMPI_Initの前に行う場合、以下のような実装が可能です。

#include <mpi.h>
#include <openacc.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
int main(int argc, char *argv[]) {
    const char* local_rank_env = getenv("OMPI_COMM_WORLD_LOCAL_RANK");
    assert(local_rank_env != NULL);
    const char* local_size_env = getenv("OMPI_COMM_WORLD_LOCAL_SIZE");
    assert(local_size_env != NULL);
    // gpu_id = [0, 1, 2, ..., ngpu_per_node - 1]
    const int gpu_id = atoi(local_rank_env);
    const int ngpu_per_node = atoi(local_size_env);
    acc_set_device_num(gpu_id, acc_device_nvidia);
    MPI_Init(&argc, &argv);
    int rank, nprocs;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
    printf("Rank %d/%d = GPU %d/%d\n", rank, nprocs, gpu_id, ngpu_per_node);
    MPI_Finalize();
}

OpenACC/MPI実装

上記のOpenACCとMPIを使ったマルチGPUの利用コードを使って、saxpyを実装してみます。
ここではMPIのsaxpy実装をOpenACC化し、OpenACC/MPI実装を作成します。

ベクトル計算をMPIで並列化するには、ベクトルをMPIプロセス数で分割し、分割したデータを各プロセスに分散処理します。

MPI実装を以下のコードに示します。最後の出力ではMPIランクの順番で計算したデータを表示しています。
MPI_Barrierで同期をしない場合、出力がランクごとにバラバラに行われるため、出力のソートの代わりに利用しています。

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <vector>
int main(int argc, char *argv[]) {
    MPI_Init(&argc, &argv);
    int rank, nprocs;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
    constexpr std::size_t N = 16;
    float alpha = 0.2f;
    const std::size_t N_per_process = (N + nprocs - 1) / nprocs;
    assert(N % nprocs == 0);
    std::vector<float> x(N_per_process), y(N_per_process);
    for (std::size_t i = 0 ; i < N_per_process ; ++i) {
        x[i] = float(rank * N_per_process + i);
        y[i] = 0.0f;
    }
    for (std::size_t i = 0 ; i < N_per_process ; ++i)
        y[i] = alpha * x[i] + y[i];
    // print sequential
    for (int j = 0 ; j < nprocs ; ++j) {
        if (j == rank)
            for (std::size_t i = 0 ; i < N_per_process ; ++i)
                printf("y[%2d] = %f (rank = %d)\n", rank * N_per_process + i, y[i], rank);
        MPI_Barrier(MPI_COMM_WORLD);
    }
    MPI_Finalize();
}

これを4プロセスで実行すると、以下の出力が得られます。
標準出力は一般的にバッファリングされているため、計算機によっては出力順序が変わることがあります。

$ mpicxx -o ./xtest mpi_axpy.cc
$ mpirun -n 4 --bind-to none ./xtest
y[ 0] = 0.000000 (rank = 0)
y[ 1] = 0.200000 (rank = 0)
y[ 2] = 0.400000 (rank = 0)
y[ 3] = 0.600000 (rank = 0)
y[ 4] = 0.800000 (rank = 1)
y[ 5] = 1.000000 (rank = 1)
y[ 6] = 1.200000 (rank = 1)
y[ 7] = 1.400000 (rank = 1)
y[ 8] = 1.600000 (rank = 2)
y[ 9] = 1.800000 (rank = 2)
y[10] = 2.000000 (rank = 2)
y[11] = 2.200000 (rank = 2)
y[12] = 2.400000 (rank = 3)
y[13] = 2.600000 (rank = 3)
y[14] = 2.800000 (rank = 3)
y[15] = 3.000000 (rank = 3)

このようなMPI化されたプログラムがあったとき、GPU化するのは非常に簡単です。
先に示したようにMPIの各ランクが使うGPUを決定した後、CPU-GPU間のデータ転送と計算カーネルのOpenACC化だけです。

以下がOpenACC/MPI実装で、前述のMPI実装と大きく変わらないことがわかると思います。
MPI実装と同じように実行してみると、同じ出力が得られることを確認してみてください。

#include <mpi.h>
#include <openacc.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <vector>
int main(int argc, char *argv[]) {
    MPI_Init(&argc, &argv);
    int rank, nprocs;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
    const int ngpu_per_node = acc_get_num_devices(acc_device_nvidia);
    assert(ngpu_per_node > 0);
    // gpu_id = [0, 1, 2, ..., ngpu_per_node - 1]
    const int gpu_id = rank % ngpu_per_node;
    acc_set_device_num(gpu_id, acc_device_nvidia);
    constexpr std::size_t N = 16;
    float alpha = 0.2f;
    const std::size_t N_per_process = (N + nprocs - 1) / nprocs;
    assert(N % nprocs == 0);
    std::vector<float> x(N_per_process), y(N_per_process);
    float *px = x.data(), *py = y.data();
    for (std::size_t i = 0 ; i < N_per_process ; ++i) {
        px[i] = float(rank * N_per_process + i);
        py[i] = 0.0f;
    }
#pragma acc data copyin(px[0:N_per_process]) copy(py[0:N_per_process])
#pragma acc kernels
    for (std::size_t i = 0 ; i < N_per_process ; ++i)
        py[i] = alpha * px[i] + py[i];
    // print sequential
    for (int j = 0 ; j < nprocs ; ++j) {
        if (j == rank)
            for (std::size_t i = 0 ; i < N_per_process ; ++i)
                printf("y[%2d] = %f (rank = %d)\n", rank * N_per_process + i, y[i], rank);
        MPI_Barrier(MPI_COMM_WORLD);
    }
    MPI_Finalize();
}

またdiffを取って確認してみると、約13行の変更だけでOpenACC/MPI化ができていることがわかります。
MPIによる分散並列はマルチGPUの場合でもそのまま利用できるので、非MPIコードをOpenACC化するのと殆ど変わらないことがわかるかと思います。

@@ -1,4 +1,5 @@
 #include <mpi.h>
+#include <openacc.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <assert.h>
@@ -11,6 +12,13 @@
     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
     MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+    const int ngpu_per_node = acc_get_num_devices(acc_device_nvidia);
+    assert(ngpu_per_node > 0);
+
+    // gpu_id = [0, 1, 2, ..., ngpu_per_node - 1]
+    const int gpu_id = rank % ngpu_per_node;
+    acc_set_device_num(gpu_id, acc_device_nvidia);
+
     constexpr std::size_t N = 16;
     float alpha = 0.2f;
@@ -18,14 +26,17 @@
     assert(N % nprocs == 0);
     std::vector<float> x(N_per_process), y(N_per_process);
+    float *px = x.data(), *py = y.data();
     for (std::size_t i = 0 ; i < N_per_process ; ++i) {
-        x[i] = float(rank * N_per_process + i);
-        y[i] = 0.0f;
+        px[i] = float(rank * N_per_process + i);
+        py[i] = 0.0f;
     }
-    for (std::size_t i = 0 ; i < N_per_process ; ++i)
-        y[i] = alpha * x[i] + y[i];
+#pragma acc data copyin(px[0:N_per_process]) copy(py[0:N_per_process])
+#pragma acc kernels
+    for (std::size_t i = 0 ; i < N_per_process ; ++i)
+        py[i] = alpha * px[i] + py[i];
     // print sequential
     for (int j = 0 ; j < nprocs ; ++j) {

まとめ

実装したMPIプログラムの計算すべてがOpenACC化されていれば、IO処理を除いてCPU-GPU間のメモリコピーは計算全体の最初と最後だけになると期待されます。
したがって、理想的なOpenACC/MPIプログラムは以下のコード構造を持つはずです。

• MPIの初期化
• 計算するデータ初期値を構築
• CPUで構築したデータをGPUへ転送
• 計算処理（時間発展など何かしらの反復処理）
• GPUからCPUへ計算結果の転送
• データの書き出しなどクリーンアップ処理

この記事で示したコードはシンプルな計算のため、ありがたみがわかりにくいかもしれません。
手元にあるプログラムはどのぐらいのコード量になっていて、それをOpenACCやMPIに対応させるとどの程度の労力がかかるかを考えるとわかりやすいと思います。

昨今はGPUのメモリ容量も増え、エンタープライズ向けではGPUあたり40–80 GiBのメモリがあるため、容量不足で実行できないことは減ったかと思います。
しかし計算量が莫大になり1 GPUでは現実的な時間で完了しない、それでもメモリが足りない、といった問題ではマルチGPU化も1つの選択肢に入ってくるでしょう。
複数台のGPUが接続された環境がある場合、OpenACCのマルチGPU化のほうがより簡単です。
一方、利用できる計算機1台が持つGPUでは足りない、将来的に複数の計算機のGPUをまとめて使う予定がある、などの理由であればOpenACC/MPIは唯一の選択肢といえます（一般的な方法の中では）。

CUDA/MPIよりは低コストですが、OpenACC/MPI化も相当のコストがかかりますので、利用する場合は費用対効果をよく検討することをおすすめします。