C++でOpenMPを使って並列化したい。

Question

Visual Studio2019でC++を使ってプログラムを書いています。

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include<time.h>

int main() {
	int i, j;
	double k = 0;
	clock_t start = clock();
#pragma omp parallel for
	for (i = 0; i <1000000 ; i++) {
		for (j = 0; j < 100; j++) {
#pragma omp atomic
		k += sqrt(i)-j;
		}
	}
	clock_t end = clock();

	const double time = static_cast<double>(end - start) / CLOCKS_PER_SEC * 1000.0;
	printf("time %lf[ms]\n", time);
	printf(" k = %lf,\n", k);
}

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include<time.h>

int main() {
	int i, j;
	double k = 0;
	clock_t start = clock();
	for (i = 0; i <1000000 ; i++) {
		for (j = 0; j < 100; j++) {
		k += sqrt(i)-j;
		}
	}
	clock_t end = clock();

	const double time = static_cast<double>(end - start) / CLOCKS_PER_SEC * 1000.0;
	printf("time %lf[ms]\n", time);
	printf(" k = %lf,\n", k);
}

mp.cppが並列化したプログラムで、notmp.cppが同じ内容で非並列化のプログラムです。
実際に実行してみた出力結果は、

time 17765.000000[ms]
k = 61624080556.706154

time 1515.000000[ms]
k = 61716616645.940987

となって、並列化しない場合のほうが10倍近く高速であるという結果になりました。
また、計算結果についても、notmp.cppでは何度やっても全く同じ結果が出力されるのに対して、
mp.cppではプログラムを動かすたびに違う結果が出てきます。
（誤差ではある。例えば61656649114.643929）
以上、余計に時間がかかる。計算結果が正確ではない。の二つの問題があり、この原因と解決策を知りたくて質問させていただきました。どうかよろしくお願いします。

Comments

[tmp]

結果のｋがかなり大きくなっているので最後の方が情報落ちが多くなりますよね。mp.cppでは、足す順番がかわるので、その誤差も変わると考えられます。

Answer 1

OpenMPもしばらく触ってない上に、コードも動かしてないのであくまで想像の上での回答で申し訳ないです。

CのOpenMPのpararell指示句で暗黙的にprivate変数になるのは並列化されたループの制御変数だけだったと思います(fortranはもうちょっと気が利いた気がしますけど)。
そのため内側のループのjはprivate変数となっておらず、各スレッドが共有する変数となってしまい、実行するたびに結果が異なってしまう原因になっていると考えられます。
またatomic指示句には若干のオーバヘッドが伴います。そのため並列数やループ数にもよりますが、オーバヘッドが無視できない状態になっているのではないかと考えられます。atomicやcriticalを安易に使用すると性能が出ない、なんていうことはよくあることなので…。こういった計算結果を足し合わせて最終的に1つの計算結果を出したい、という場合にはreduction指示句を使用する場合が多いのではないかと思います。以下のページが参考になるのではないでしょうか。

もくもくブログ　OpenMPによるfor文の並列化

Answer 2

余計に時間がかかる。計算結果が正確ではない。

Munosuke222さんの回答内容そのままですが、下記のように変更したら、notmp.cppよりも速く、実行する度に計算結果が変わらなくなりました
ただし、notmp.cppの計算結果と同じではありませんでした (違いは1.0未満ですが)

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include<time.h>

int main() {
	int i, j;
	double k = 0;
	clock_t start = clock();
//#pragma omp parallel for
//#pragma omp parallel for private(i, j)
#pragma omp parallel for private(i, j) reduction(+:k)
	for (i = 0; i <1000000 ; i++) {
		for (j = 0; j < 100; j++) {
//#pragma omp atomic
		k += sqrt(i)-j;
		}
	}
	clock_t end = clock();

	const double time = static_cast<double>(end - start) / CLOCKS_PER_SEC * 1000.0;
	printf("time %lf[ms]\n", time);
	printf(" k = %lf,\n", k);
}

　
Visual Studio 2017で、下記のコマンドでコンパイルしました

cl /O2 /arch:AVX /fp:fast /EHsc /MT /openmp /wd 4819 mp.cpp

Comments

[jbpb0]

当方のPCでの結果は、下記の通り

[mp.cppの場合]
time 13.000000[ms]
k = 61716616645.882156,

[notmp.cppの場合]
time 55.000000[ms]
k = 61716616645.764549,

コンパイルコマンドに「/fp:fast」を付けなかった場合の結果は下記の通りで、notmp.cppの結果の「k」は、質問者さんの結果と同じでした

[mp.cppの場合]
time 17.000000[ms]
k = 61716616645.887421,

[notmp.cppの場合]
time 118.000000[ms]
k = 61716616645.940987,

【追記】
コンパイルコマンドに「/O2 /arch:AVX /fp:fast」を付けなかった場合の結果は下記の通り

[mp.cppの場合]
time 150.000000[ms]
k = 61716616645.887421,

[notmp.cppの場合]
time 1124.000000[ms]
k = 61716616645.940987,

[Munosuke222]

これくらいの計算だと並列化とベクトル化で随分と速くなるんですねえ…