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)で定義した演算を使う
追記:
- 要素
ans[i]自体に重複書き込みする可能性がなければ、reductionの記述は不要- この例では不要だった
-
#pragma omp parallel forと1行で書くこともできる
書き込み変数が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