5. OpenMP

5.1 OpenMPとは

OpenMP[3][4]とは、共有メモリー環境で並列計算を行うプログラミング技術であり、 複数のコアを持ったCPUで、 同時に複数のスレッドを実行して処理を分割し計算時間を短縮するものです。 1台のコンピュータに複数のCPUが搭載されていればさらにその個数だけ速くなります。
OpenMPは次節で述べるマルチスレッドの一種ですが、 プログラミング方法が異なりますので、ここでは特に"OpenMP"と呼びます。

5.2 OpenMPプログラミング

OpenMPでは並列処理を行うfor文のすぐ上に"#pragma omp"で始まる指示文(ディレクティブ)を挿入し、 コンパイルオプションでOpenMPを使用していることを知らせます。
コンパイルオプションがないときは指示文はコメントとみなされ 1コアのみを使用する逐次プログラムができます。
OpenMPは元のソースコードに一切手を加えることなく並列化できますので、 OpenMPで並列化できるアルゴリズムのときは最も簡単にプログラミングできます。 また、並列化したい場所だけを順々に並列化できますのでプログラミングの生産性が高くなります。
なお、OpenMPはOSに依存しない規格ですので、同じソースコードがWindowsとLinuxで動きます。

5.3 OpenMPプログラミング例

リスト5-1にベクトルの和または内積をOpenMPで並列計算するプログラムを示します。


リスト5-1 OpenMPプログラム(omp.c)
     1	/*
     2	add vector, dot product (OpenMP)
     3	
     4	VC++ : cl /Ox /openmp omp.c
     5	gcc  : gcc -O3 -std=c99 -fopenmp omp.c -o omp
     6	
     7	Usage :
     8	> omp {1|2} <num> <loop> <thread>
     9	*/
    10	
    11	#include <stdlib.h>
    12	#include <stdio.h>
    13	#include <time.h>
    14	
    15	#ifdef _OPENMP
    16	#include <omp.h>
    17	#endif
    18	
    19	#ifndef _WIN32
    20	#include <sys/time.h>
    21	#endif
    22	
    23	static void vadd(int n, const float a[], const float b[], float c[])
    24	{
    25		int    i;
    26	#ifdef _OPENMP
    27	#pragma omp parallel for
    28	#endif
    29		for (i = 0; i < n; i++) {
    30			c[i] = a[i] + b[i];
    31		}
    32	}
    33	
    34	static double sdot(int n, const float a[], const float b[])
    35	{
    36		double sum = 0;
    37		int    i;
    38	#ifdef _OPENMP
    39	#pragma omp parallel for reduction(+:sum)
    40	#endif
    41		for (i = 0; i < n; i++) {
    42			sum += a[i] * b[i];
    43		}
    44	
    45		return sum;
    46	}
    47	
    48	int main(int argc, char **argv)
    49	{
    50		int    mode = 1;
    51		int    nthread = 1;
    52		int    n = 1000;
    53		int    nloop = 1000;
    54	#ifdef _WIN32
    55		clock_t t0, t1;
    56	#else
    57		struct timeval t0, t1;
    58	#endif
    59	
    60		// arguments
    61		if (argc >= 5) {
    62			mode    = atoi(argv[1]);
    63			n       = atoi(argv[2]);
    64			nloop   = atoi(argv[3]);
    65			nthread = atoi(argv[4]);
    66		}
    67	
    68	#ifdef _OPENMP
    69		// thread
    70		omp_set_num_threads(nthread);
    71	#endif
    72	
    73		// alloc
    74		size_t size = n * sizeof(float);
    75		float *a = (float *)malloc(size);
    76		float *b = (float *)malloc(size);
    77		float *c = (float *)malloc(size);
    78	
    79		// setup problem
    80		for (int i = 0; i < n; i++) {
    81			a[i] = i + 1.0f;
    82			b[i] = i + 1.0f;
    83		}
    84	
    85		// timer
    86	#ifdef _WIN32
    87		t0 = clock();
    88	#else
    89		gettimeofday(&t0, NULL);
    90	#endif
    91	
    92		// calculation
    93		double sum = 0;
    94		if (mode == 1) {
    95			// add vector
    96			for (int loop = 0; loop < nloop; loop++) {
    97				vadd(n, a, b, c);
    98			}
    99			for (int i = 0; i < n; i++) {
   100				sum += c[i];
   101			}
   102		}
   103		else {
   104			// dot product
   105			for (int loop = 0; loop < nloop; loop++) {
   106				sum += sdot(n, a, b);
   107			}
   108			sum /= nloop;
   109		}
   110	
   111		// timer
   112		double cpu = 0;
   113	#ifdef _WIN32
   114		t1 = clock();
   115		cpu = (double)(t1 - t0) / CLOCKS_PER_SEC;
   116	#else
   117		gettimeofday(&t1, NULL);
   118		cpu = (t1.tv_sec - t0.tv_sec) + 1e-6 * (t1.tv_usec - t0.tv_usec);
   119	#endif
   120	
   121		// output
   122		double exact = (mode == 1)
   123			? (double)n * (n + 1)
   124			: (double)n * (n + 1) * (2 * n + 1) / 6.0;
   125		printf("n=%d nloop=%d nthread=%d sum=%.6e(%.6e) cpu[sec]=%.3f\n",
   126			n, nloop, nthread, sum, exact, cpu);
   127	
   128		// free
   129		free(a);
   130		free(b);
   131		free(c);
   132	
   133		return 0;
   134	}

OpenMP指示文
コンパイルオプション/openmpをつけるとマクロ"_OPENMP"が有効になります。 27行目と39行目でfor文を並列計算することを指示しています。 なお、39行目ではベクトルの内積を並列計算するには、 すべてのスレッドが総和変数にアクセスしますので以下の"reduction"指示文が必要になります。
reduction (operatorlist)
ここで operator は加算("+")、乗算("*")などを表し、 list はreduction演算される変数名を表します。

OpenMP関数
70行目のomp_set_num_threads関数は使用するスレッド数を指定する関数です。 本関数を呼ばないときは論理コア数が指定されます。
OpenMP関数を使用するには16行目のinclude文が必要になります。
OpenMPにはこの他にいくつかの関数があります。 使用頻度は高くありませんが、詳しくは[3][4]を参考にして下さい。

計算時間の計測
gccではclock関数は全スレッドの和になりますので、 リスト5-1では計算時間を計測するためにgettimeofday関数を使用しています。

コンパイル・リンク方法
VC++ではコンパイルオプションに"/openmp"が必要です。
> cl /Ox /openmp omp.c
gccではコンパイルオプションとリンクオプションに"-fopenmp"が必要です。
> gcc -O3 -std=c99 -fopenmp omp.c -o omp

プログラムの実行方法
プログラムの実行方法は以下の通りです。
> omp {1|2} 配列の大きさ 繰り返し回数 スレッド数
例えば以下のようになります。
> omp 1 1000000 1000 4 (ベクトルの和を計算するとき)
> omp 2 1000000 1000 4 (ベクトルの内積を計算するとき)
繰り返し回数は計算時間の測定誤差を小さくするためです。
なお、OpenMP対応プログラムの実行にはvcomp140.dllファイルが必要です。 (Microsoft Visual Studio 2017の場合)

5.4 ベクトル和の計算時間

表5-1にベクトル和の計算時間を示します。
配列の大きさ(=N)と繰り返し回数(=L)の積は一定(=1010)です。従って全体の演算量は同じです。
表の下になるほど使用メモリーが小さくなりキャッシュに乗りやすくなります。
No.3のとき最も並列計算の速度比が高くなっています。

表5-1 ベクトル和の計算時間(OpenMP, Windows, ()内は1スレッドとの速度比)
No.配列の大きさN繰り返し回数L 1スレッド 2スレッド 4スレッド 8スレッド
1 10,000,000 1,00011.06秒(1.0)11.00秒(1.01)11.12秒(0.99)11.22秒(0.99)
2 1,000,000 10,000 7.22秒(1.0) 6.56秒(1.10) 6.04秒(1.20) 5.89秒(1.23)
3 100,000 100,000 4.38秒(1.0) 2.32秒(1.89) 1.53秒(2.86) 1.23秒(3.56)
4 10,000 1,000,000 4.01秒(1.0) 2.52秒(1.59) 2.30秒(1.74) 2.11秒(1.90)

5.5 ベクトル内積の計算時間

表5-2にベクトル内積の計算時間を示します。
配列の大きさ(=N)と繰り返し回数(=L)の積は一定(=1010)です。従って全体の演算量は同じです。
表の下になるほど使用メモリーが小さくなりキャッシュに乗りやすくなります。
No.2-3のとき最も並列計算の速度比が高くなっています。

表5-2 ベクトル内積の計算時間(OpenMP, Windows, ()内は1スレッドとの速度比)
No.配列の大きさN繰り返し回数L 1スレッド 2スレッド 4スレッド 8スレッド
1 10,000,000 1,0008.71秒(1.0)6.09秒(1.43)5.82秒(1.50)5.94秒(1.47)
2 1,000,000 10,0008.03秒(1.0)4.11秒(1.95)2.44秒(3.29)1.54秒(5.21)
3 100,000 100,0007.80秒(1.0)4.08秒(1.91)2.53秒(3.08)1.59秒(4.91)
4 10,000 1,000,0007.80秒(1.0)4.61秒(1.69)3.15秒(2.48)2.82秒(2.77)