Edited at

C++並列処理: openmp, thread, Eigen


C++並列処理: openmp, thread, Eigen

毎回忘れてしまうC++の並列処理のメモ.

よく使うやつだけ書きます.


Openmp

おそらく最も単純で誰でもすぐに使える並列化の方法.

コンパイル時に, -fopenmpオプションをつける必要がある.


forを並列化する


書き込み変数がスカラーの場合

ベクトルの内積を求める例で示す.

double dot(const vector<double> &v1, const vector<double> &v2) {

double ans = 0;
#pragma omp parallel
{
#pragma omp for reduction(+:ans)
for (size_t i = 0; i < v1.size(); ++i){
ans += v1[i]*v2[i];
}
}
return ans;
}



  • #pragma omp parallelで並列化対象を囲む


  • #pragma omp forでforを並列化する

  • CPUが同時に書き換える恐れのある変数ansは, reductionとしてスレッドごとに別々に扱うようにする

  • 最終的に, スレッドごとのans+でまとめればよいため, reduction(+:ans)と書く


書き込み変数がVectorの場合

同時に書き換える恐れのある変数がvectorの場合は少し面倒になる.

ベクトル同士の要素ごとの積を例に示す.

#pragma omp declare reduction(vec_double_plus : std::vector<double> : \

std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), std::plus<double>())) \
initializer(omp_priv = omp_orig)
vector<double> cwiseProduct(const vector<double> &v1, const vector<double> &v2) {
vector<double> ans(v1.size(),0);
#pragma omp parallel
{
#pragma omp for reduction(vec_double_plus :ans)
for (size_t i=0;i<v1.size();++i){
ans[i]=v1[i]*v2[i];
}
}
return ans;
}



  • #pragma omp declare reductionでreductionで使う演算について定義する


  • vec_double_plus : std::vector<double>で, vector<double>を引数にとる演算vec_double_plusを定義する(2つのベクトルの要素ごとの和)



    • std::transformには#include <algorithm>が必要




  • #pragma omp for reduction(vec_double_plus :ans)で定義した演算を使う


書き込み変数がMatrixの場合

EigenのMatrixXdを例に書く.

MatrixXd A=MatrixXd::Zero(100,100);

#pragma omp declare reduction(+ : Eigen::MatrixXd : omp_out=omp_out+omp_in) initializer(omp_priv = omp_orig)
#pragma omp parallel for reduction(+:A)
for (int i=0;i<10;++i){
A+=MatrixXd::Random(100,100);
}


  • Eigenはoperator+が定義されているが, #pragma omp declareしないといけない


  • std::vectorよりもこちらのほうがわかりやすいのは, operator+がEigenで定義されているため


    • そのため, std::vectorの方も演算子+を定義してからの方がわかりやすいかもしれない




  • #pragma omp parallel for#pragma omp parallel#pragma omp forの省略表現


タスクの並列化

複数のタスクを同時に実行したいときに有効.

#pragma omp parallel

{
#pragma omp sections
{
#pragma omp section
{
dot(v1,v2);
}
#pragma omp section
{
cwiseProduct(v1,v2);
}
}
}



  • #pragma omp parallelで並列化範囲を宣言


  • #pragma omp sectionsでセクションごとで並列化することを宣言


  • #pragma omp sectionで1つのセクションを宣言する


Thread

コンパイル時に-pthreadが必要.

#include<thread>が必要.


戻り値なしの場合

以下の例では, 関数addRefの実行N回を別プロセスで行う.

結果を受け取るには参照渡しをするしかない.

void addRef(int x, int y,int &ans){

ans = x + y;
}
int main() {
const int N = 5;
const int x=1, y=2;

vector<thread> threads;
vector<int> result(N);

for (int i=0;i<N;++i)
threads.push_back(thread(addRef, x, y, ref(result[i])));
for (thread &t: threads)
t.join();
for (int r : result)
cout<<r<<endl;
return 0;
}


  • 最初のforでスレッド作成: thread(関数名, 引数1, 引数2, ...)



    • refを使わないと参照渡しできない



  • 次のforでスレッド終了待ち: join()

  • 最後のforで結果出力


戻り値ありの場合

ラムダ式で書けば受け取れる.

int addReturn(int x, int y){

return x + y;
}
int main() {
const int N = 5;
const int x=1, y=2;

vector<thread> threads;
vector<int> result(N);

for (int i=0;i<N;++i)
threads.push_back(thread([i, &result]{ result[i] = addReturn(x, y); }));
for (thread &t: threads)
t.join();
for (int r:result)
cout<<r<<endl;
return 0;
}


  • 最初のforをlambda式で覆っただけ


  • asyncを使ってもいい


Eigen

Eigenは, 行列を扱うライブラリ.

ヘッダーをインクルードするだけで使える: (例)g++ -I./Eigen test.cpp#include <Eigen/Dense>

クイックリファレンス: AsciiQuickReference.txt

Eigenは-fopenmpオプションをつけるだけで自動的に並列化される.

公式ドキュメントMake Eigen run in parallelを見ると, 以下の機能は並列化に対応しているようだ.


  • general dense matrix - matrix products

  • PartialPivLU

  • row-major-sparse * dense vector/matrix products

  • ConjugateGradient with Lower|Upper as the UpLo template parameter.

  • BiCGSTAB with a row-major sparse matrix format.

  • LeastSquaresConjugateGradient


参考文献