Physics-Station 2

旧 http://physics-station.blogspot.jp/ から当はてなブログに移行しました。間違ってるところがあればコメントください。記述の正確性は保証しません。

OpenMPでfor文を高速化する

openmpでfor文を高速化してみよう。OpenMPを有効にするには、Visual Studioのプロジェクトのプロパティページを開いて、C/C++の言語のOpenMPのサポートをはい(/openmp)にする必要がある。

例として1から50000までの数が素数(prime)かどうかを調べてみよう。素数の調べ方は2以上候補値未満の数で割った余りを調べる愚直な方法で実装した。

#include <omp.h>
#include <vector>
#include <iostream>
#include <chrono>
#include <algorithm>

int main()
{
    auto *mylock = new omp_lock_t;
    omp_init_lock(mylock); //lockを初期化

    std::vector<int> result; //素数の結果

    auto start = std::chrono::system_clock::now(); //開始時間

#pragma omp parallel for num_threads(12)
    for (int i = 0; i < 100; i++) {
        std::vector<int> list; //素数の中間ファイル
        
        for (int j = 0; j < 5000; j++) {
            int candidate = i + j * 100 + 1;
            if (candidate < 2) { continue; }

            bool prime_flag = true;
            for (int k = 2; k < candidate; k++) {
                if ((candidate%k) == 0) { //2以上候補値未満の数で割った余り
                    prime_flag = false; //素数ではない
                    break;
                }
            }
            if (prime_flag) {
                list.emplace_back(candidate);
            }
        }
        omp_set_lock(mylock); //std::vector<T>の変更操作をするためlock
        for (auto p = list.begin(); p != list.end(); p++) {
            result.emplace_back(*p);
        }
        omp_unset_lock(mylock); //lockを解除
    }
    auto end = std::chrono::system_clock::now(); //終了時間

    omp_destroy_lock(mylock); //lock破棄

    std::cout << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << " [msec]" << std::endl;

    std::cout << "max prime is " << *std::max_element(result.begin(), result.end()) << std::endl; //調査した中で最大の素数を出力
    ::system("pause");
}

当環境で 37秒かかったものが12スレッド並列で行うと6.5秒に短縮された。オーバーヘッドがあるため、単純に1/12にはならない。 std::vector<T>の変更操作はスレッドセーフではないため、omp_lock_tを使ってlockとlock解除を行っている。

素数の調べ方は最適化しておらず、世の中にはもっと高速な方法が沢山ある。一般的に、openmpで高速化するのはアルゴリズムを改善した後のほうが良いだろう。

スレッド数を指定しない場合は以下のように記述することも可能。この場合は全スレッドを使って計算が行われる。

#pragma omp parallel for