この記事は IPFactory OB Advent Calendar 2023 15日目の記事です。

こんにちは、ふたばとです。

富岳のようないわゆるスパコンでパソカタをしていくには、並列計算の技法に関する知識が必要になってきます。本稿では、並列計算の概要と OpenMP を使ってマルチスレッドでの並列計算を見ていきます。

手っ取り早く並列計算をしたい場合には OpenMP でシュッとするのはいいのですが、大規模に並列させるには MPI(Multi Process Interface) を使うのが基本になります。今回は OpenMP のマルチスレッドの話にはなりますが、近い機会にアクセラレータを使う OpenACC や MPI についての記事も書いてみようと思います。

なお、本稿では intel コンパイラ2022.3.1 の環境で実行させています。

並列計算について

並列計算とは、同時に計算資源を利用して計算を行う手法です。演算を分割して複数のコアやプロセッサに割り当てて処理させることで実現されます。並列計算の手法として、スレッド並列、プロセス並列、アクセラレータを用いた並列化があります。それぞれの手法には特有の利点や制約があり、計算機アーキクチャ、ハードウェア、アプリケーションに依存するので、計算機や解く問題の特性に合わせたプログラミングが求められます。

スレッド並列（共有メモリ型）

スレッド並列型の並列計算は、分割された演算（スレッド）がメモリ空間を共有できることが特徴になります。複数のスレッドが均一に並列処理をするため、メモリ空間は単一で、手軽に並列化の実装が可能です。

基本はマスタースレッドだけで計算させて、ある特定の演算（負荷の多いループ・行列演算など）に対してスレッドチームを組んで複数のスレッドで処理をさせます。実装には並列で計算してほしいところに並列処理を実装すればよいです。なおプロセス間の通信にはオーバーヘッドがあるため、適切なデータ同期が必要にはなります。

プロセス並列（分散メモリ型）

プロセス並列型の並列計算は、分割された演算（プロセス）がそれぞれ独立のメモリ空間を参照します。全プロセスが同一のプログラムを走らせるのですが、それぞれのプロセスで持っているデータが違うために出力も変わってきます。プロセス並列型は自身のメモリ空間上のデータのみアクセス可能なため、他のプロセスのデータへのアクセスには MPI 通信による転送が必要で、プロセス間のデータ通信を明示的に記述する必要があります。また、各 MPI プロセスは固有のランク番号を持ち、各プロセスで処理を変えていきたい場合には、ランク番号を利用した条件分岐を記述する必要があります。スレッド並列型の並列計算手法に比べて、プロセス並列型はやや実装の難易度は高くなります。

なお、富岳をはじめとした大規模なスパコンは大抵ハイブリッドな構成になっています。大規模な CPU を用意するにはまず分散メモリ型で計算機を構成する必要があり、さらに一つの計算機の中で共有メモリ型の計算機を構成していきます。よって共有メモリ型の並列計算と分散メモリ型の並列計算の両方を知っておくとよいと考えられます。

アクセラレータを用いた並列化

まず、アクセラレータは CPU の処理を一部代替して全体の処理の高速化を図る装置です。 CPU（ホスト）からPCI バスを通じてアクセラレータ（デバイス）のメモリにアクセスします。デバイスの計算で必要なデータはホストから転送し、デバイス側の処理が終わればホストに戻す処理を記述する必要があり、両者のメモリ空間は異なるため、うまく高速化できるような並列処理の実装は慣れが必要で、個人的には簡単ではないなと感じました。

OpenMPでマルチスレッドな並列計算をシュッと実装する

OpenMP（Open Multi-Processing）は、共有メモリ型並列計算機環境環境においてマルチスレッド並列処理を行うための API です。 OpenMP はC、C++、Fortranなどのプログラミング言語で利用可能です。

OpenMP はディレクティブベースのプログラミンング手法を採用しています。つまりディレクティブ（プリプロセッサによって認識される特定の命令）を使用して並列化の指示を行います。実装する人間は、並列で計算してほしいところに並列処理を実装すればよいため、元のコードを大きく変更することなく並列化処理を導入できます。

OpenMP の基本

おまじないのロードは #include <omp.h> です。

#include <stdio.h>
#include <omp.h>

int main()
{
    int i;

    // OpenMP を有効にするディレクティブ
    #pragma omp parallel
    {
        // 並列実行領域
        #pragma omp for
        for(i=0; i<1000; i++)
        {
            // ループが並列処理される
        }
    }

    return 0;
}

parallel で囲われている部分が並列実行領域となります。

最大スレッド数の取得に omp_get_num_threads()、自スレッド番号の取得に omp_get_thread_num() が使えます。並列計算で利用するスレッド数は環境変数 OMP_NUM_THREADS として設定します。コンパイルには -qopenmp オプションを付与して実行します。-qopenmp オプションを付与しない場合、指示文はコメントとして扱われます。

$ icx -qopenmp hoge.c

データスコープ属性

マルチスレッド型の並列計算はメモリを共有しています。メモリを共有しているため、各スレッドは共有メモリ上の変数にアクセスすることができます。故にすべての変数をグローバルにアクセル可能であると、スレッド同士がお互いの処理を邪魔してしまう可能性があります。このため、複数のスレッドによる同時のアクセスや変更が発生しないようにするために、変数にデータスコープ属性を付与してプログラム内での変数のスコープとアクセス範囲を制御することが重要になります。

概念的には全く難しくないと思いますが、スコープを意識してプログラミングをしないとバグを引き起こす原因となるため、注意が必要です。

shared 変数

shared 変数はすべてのスレッドで同一のメモリを参照します。複数のスレッドが同じデータを更新する必要がある場合、shared 変数を使用することが適しています。また、大規模なデータセットに対して同じ処理を複数のスレッドで行う場合も shared 変数として扱うこともあるでしょう。

ただし、shared 変数を使用する際には注意が必要で、競合状態（Race Condition）を引き起こす可能性があるためデータの整合性の確保が求められます。競合状態は、複数のスレッドが同じメモリ位置に対して同時に書き込む場合や、一方が書き込みを行っている最中に他方が読み込みを行う場合などに発生します。競合状態が発生すると、プログラムの動作が不安定になり、予測できない結果が生じる可能性があります。

private 変数

private 変数は各スレッドが異なるメモリを参照します。他のスレッドが上書きしてしまうのが怖いループの index や一時変数は private 変数として指示します。 private 変数の値はスレッドごとに異なる値になっていることが期待されるため、その特性を理解して利用する必要があるでしょう。

暗黙のデータスコープ属性

データスコープ属性が指定されていない変数、並列処理の領域外で定義されている変数は shared 変数です。

ループの index や並列処理の領域内で定義されている変数は private 変数になります。ただし、多重ループの場合、内側のループ（以下の例でいう j, k のループ）は shared 変数扱いのため、private 変数であることを明示する必要があります。

#pragma omp parallel
{
    #pragma omp for private(j, k)
    for (i=0; i<1000; i++) {
        for (j=0; j<1000; j++) {
            for (k=0; k<1000; k++) {
                a[i][j][k] = b[i] * c[j] * d[k];
            }
        }
    }
}

データスコープのデフォルト動作を制約する

通常、OpenMPディレクティブ内で変数を使用する際、その変数が shared 変数であるか、private 変数であるかを指定しますが、default(none) を指定すると、全ての変数が明示的に指定されない限り、コンパイラがエラーを発生させるようになります。

#pragma omp parallel default(none) shared(x) private(y) 
{
    // x is shared, y is private
}

default(none) にして各変数を明示的に shared, private と振っていく実装が好ましいです。 shared, private の振り分けを忘れてしまっていても、コンパイラが怒ってくれるようになります。

並列できるもの、並列できないもの

どのようなループが並列可能であるか、並列化に注意が必要なループには何があるかを知っておく必要があります。並列可能であるかどうかはコンパイラは教えてくれない場合があるため、基本的に人間が考える必要があります。不適切な実装によって、OpenMP は無理やり演算することで誤った演算結果を提示するでしょう。 OpenMP を使ったマルチスレッドな並列計算を実装する場合には、parallel を用いた単純なループを並列化させることが多くなると考えらえます。

ループ中に中断処理がある場合

for あるいは while ループの中に break や throw の中断処理を入れることは OpenMP が禁止しています。ループの中断を実現したい場合には、自分でフラグ管理をする必要があります。

ループ中に依存関係がある場合

以下のコードを例に、スレッドを分けて考えてみます。

for(i=0; i<1000; i++)
    a[i] = a[i-1] + b[i];

スレッド1

for(i=0; i<500; i++)
    a[i] = a[i-1] + b[i];

スレッド2

for(i=500; i<1000; i++)
    a[i] = a[i-1] + b[i];

上記のコードでは、a[i] を求めるために a[i-1] に依存していると言えます。これは a[i-1] であろうが a[i+1] であろうが同じで、前後に依存している場合には並列化はできません。

ループ中に間接参照が含まれる場合

間接参照の場合も難しいので注意が必要です。

for(i=0; i<1000; i++)
    index[i] = ......;

for(i=0; i<1000; i++)
    a[index[i]] = b[i] + b[i];

share 変数はデータの一貫性は保証しません。タイミング次第ではバグを引き起こすため、critical という排他制御ができる指示行を明示するなど注意が必要です。ただし、critical 補助指示文を利用すると高スレッド数での実行で速度が低下するようです。

ループ中に一時変数が含まれる場合

for(i=0; i<1000; i++)
{
    tmp = a[i] * b[i];
    c[i] = tmp;
}

一時変数 tmp がshared 変数の場合、異なるスレッド間で同時にアクセスされるとタイミングによって答えが違ってくることになります。ループの index を他のスレッドが持っているのは怖いので、一時変数は private 変数にして実装しましょう。

Reduction 演算

マルチスレッドな並列処理の実装には、shared 変数に対する競合状態を避けるために、適切な同期メカニズムを使用することが重要です。 shared 変数に対するスレッド間の安全なアクセスが確保のために、Reduction演算があります。 Reduction 演算は、並列処理における複数のスレッドで同時に実行される部分で、複数のスレッドが shared 変数に対して演算を行い、その結果を一つのスカラー値にまとめる操作です。

構文としては以下のように、 op が演算子 (+, - *, &, ^, |, &&, ||)、var が reduction変数です。 reduction 変数は複数の変数を指定可能です。

#pragma omp for reduction(op:var)

reduction 変数は、各スレッドの実行中には private 変数として扱われて、ループから戻ってから元の変数に集約します。

int sum = 0;
#pragma omp parallel for reduction(+:sum)
for(i=0; i<1000; i++)
    sum = sum + a[i];

#pragma omp atomic #pragma omp critical でも Reduction 演算を記述可能ですが、reduction 節を使ったほうが高速になります。ただし、reduction 節は排他的に加算が行われるものであるため、reduction 変数用の配列を用意して逐次加算していく方が高速になる場合もあるようです。

アムダールの法則

最後に、アムダールの法則について言及しておきます。

Wikipedia を眺めてもらえれば十分なのですが、アムダールの法則は並列処理における性能向上の限界を示す法則であり、並列化による効果を見積もるためのモデルです。

ja.wikipedia.org

スレッドを増やすとスピードはあがりますが、スピードの上昇率はどんどん落ちていくので、ちょうどいいスレッド数を考えてやる必要があるというものです。スレッド自体は増やせても、並列可能な割合が0.9~0.95程度でほとんど変わらないのであれば、むやみにスレッド数を増やすことは資源の使い方として、うまいとはいえないものですよねというものです。

プログラムの全体の中で直列処理（並列化できない部分）が占める割合が大きいほど、並列化による効果が制限されることを示しています。