畳み込み演算 - 視覚数理モデルシミュレーションの並列化と錯視画像の探索

function (conv̲cuda̲shared)

OpticalFlow Calculation for each image

Convolution

Data transfer GPU→CPU Data transfer

CPU→GPU

time

図

3.3.3: NVIDIA Visual Profiler

^{のタイムライン}

image1 image2

f t

S xx

u

S xt S xy S yt S yy

v

f y

f x

四則演算

1 int2 CUDATHREAD ; CUDATHREAD .x = 32; CUDATHREAD .y = 16;

2 dim3 block (( w_ex1 + CUDATHREAD .x - 1) / CUDATHREAD .x ,

3 ( h_ex1 + CUDATHREAD .y - 1) / CUDATHREAD .y ); // グリッドあたりのブロック数 =(16 ,32) 4 dim3 thread ( CUDATHREAD .x , CUDATHREAD .y ); // ブロックあたりのスレッド数 =(32 ,16) 5

6 /* それぞれの変数を格納するのに必要なメモリ領域を計算 */

7 size_t IMAGE = sizeof(float)* wh ; 8 size_t IMAGE1 = sizeof(float)* wh_ex1 ; 9 size_t IMAGE2 = sizeof(float)* wh_ex2 ; 10 size_t FILTER1 = sizeof(float)* ff1 ; 11 size_t FILTER2 = sizeof(float)* ff2 ;

12 // IMAGE1 , IMAGE2 dx ,dy , gs fx ,fy , ft fxfx ,... fxfx_ ,... rf 13 size_t MALLOC_SIZE = IMAGE *2 + IMAGE1 *2 + FILTER1 *3 + IMAGE *3 + IMAGE *10 + IMAGE2 *5 + FILTER2

+ IMAGE *13;

14 float * malloc ; cudaMalloc ((void**)& malloc , MALLOC_SIZE ); // メモリ確保 15

16 /* シェアードメモリ用の領域 */

17 int shared1 = sizeof(float)*( CUDATHREAD .x +2* side1 )*( CUDATHREAD .y +2* side1 ) + FILTER1 ; 18 int shared2 = sizeof(float)*( CUDATHREAD .x +2* side2 )*( CUDATHREAD .y +2* side2 ) + FILTER2 ; 19

20 /* 確保したメモリ領域のアドレスから，それぞれの変数の先頭部分を渡す */

21 float * img1 = & malloc [0];

22 float * img2 = & img1 [ wh ];

23 (----略- - - -) 24

25 /* C P UからG P Uへデータ転送 */

26 cudaMemcpy ( img1 , h_img1 , IMAGE , cudaMemcpyHostToDevice );

27 cudaMemcpy ( img2 , h_img2 , IMAGE , cudaMemcpyHostToDevice );

28 cudaMemcpy (dx , h_dx , FILTER1 , cudaMemcpyHostToDevice );

29 (----略- - - -) 30

31 /* 畳み込み演算 */

32 conv_kernel_shared <<< block , thread , shared1 > > >( img1 , fx , dx , w , h , f1 ); // img1 * dx → fx 33 conv_kernel_shared <<< block , thread , shared1 > > >( img1 , fy , dy , w , h , f1 ); // img1 * dy → fy 34 conv_kernel_shared <<< block , thread , shared1 > > >( img2 , ft , gs , w , h , f1 ); // img2 * gs → ft 35

36 mul_matrix_cuda <<< block , thread > > >(fx , fx , fxfx , w , h ); // fx x fx → fxfx 37 conv_kernel_shared <<< block , thread , shared2 > > >( fxfx , S_xx , rf , w , h , f2 ); // fxfx * rf → S_xx 38 mul_matrix_cuda <<< block , thread > > >(fy , fy , fyfy , w , h ); // fy x fy → fyfy 39 conv_kernel_shared <<< block , thread , shared2 > > >( fyfy , S_yy , rf , w , h , f2 ); // fxfx * rf → S_yy 40 mul_matrix_cuda <<< block , thread > > >(fx , fy , fxfy , w , h ); // fx x fy → fxfy 41 conv_kernel_shared <<< block , thread , shared2 > > >( fxfy , S_xy , rf , w , h , f2 ); // fxfy * rf → S_xy 42 mul_matrix_cuda <<< block , thread > > >(fx , ft , fxft , w , h ); // fx x ft → fxft 43 conv_kernel_shared <<< block , thread , shared2 > > >( fxft , S_xt , rf , w , h , f2 ); // fxft * rf → S_xt 44 mul_matrix_cuda <<< block , thread > > >(fy , ft , fyft , w , h ); // fy x ft → fyft 45 conv_kernel_shared <<< block , thread , shared2 > > >( fyft , S_yt , rf , w , h , f2 ); // fyft * rf → S_yt 46

47 add_matrix_cuda <<< block , thread > > >( S_xx , epsilon , w , h );

48 add_matrix_cuda <<< block , thread > > >( S_yy , epsilon , w , h );

50 mul_matrix_cuda <<< block , thread > > >( S_xx , S_yy , S_xxyy , w , h ); // S_xx x S_yy → S_xxyy 51 mul_matrix_cuda <<< block , thread > > >( S_xy , S_xy , S_xyxy , w , h ); // S_xy x S_xy → S_xyxy 52 mul_matrix_cuda <<< block , thread > > >( S_yt , S_xy , S_ytxy , w , h ); // S_yt x S_xy → S_ytxy 53 mul_matrix_cuda <<< block , thread > > >( S_xt , S_yy , S_xtyy , w , h ); // S_xt x S_yy → S_xtyy 54 mul_matrix_cuda <<< block , thread > > >( S_xt , S_xy , S_xtxy , w , h ); // S_xt x S_xy → S_xtxy 55 mul_matrix_cuda <<< block , thread > > >( S_yt , S_xx , S_ytxx , w , h ); // S_xx x S_yy → S_xxyy 56

57 sub_matrix_cuda <<< block , thread > > >( S_ytxy , S_xtyy , u , w , h ); // S_ytxy - S_xtyy → u 58 sub_matrix_cuda <<< block , thread > > >( S_xtxy , S_ytxx , v , w , h ); // S_xtxy - S_ytxx → v 59

60 (----略- - - -) 61

62 cudaDeviceSynchronize (); // スレッド同期

63 cudaMemcpy (dst , tmp_uu , IMAGE , cudaMemcpyDeviceToHost ); // G P UからC P Uへデータ転送 64 cudaFree ( malloc ); // メモリ解放

図

3.3.5:

オプティカルフロー計算の実装一部

(CUDA,

^{データ転送の削減}

)

3.3.4

マルチストリーム化

GPU

は同一時間に単一のカーネルを実行するだけではなく，計算リソースが許せば並列で複数のカーネルを実行できる．

GPU

での計算のスケジュール管理の単位をストリームと呼び，複数ストリームを利用することをマルチストリーム化と呼ぶ．マルチストリームの例を図

3.3.6

^に示す．

この機能をうまく利用するとデータ転送とカーネル実行を並行で行えたり，連続カーネル間での実行までのオーバーヘッドを減らすことが期待できる．

�me

Stream 0

cudaMemcpy

cudaMemcpy CUDA kernel CUDA kernel cudaMemcpy cudaMemcpy

Stream 0 Stream 1

cudaMemcpy

CUDA kernel CUDA kernel

cudaMemcpy

cudaMemcpy cudaMemcpy

�me

Single Stream Multi Stream

図

3.3.6:

マルチストリーム処理の例

図

3.3.4

で示したように，シミュレーションのデータ依存性が決まっている．例えば，

S

xtを求める処

理を行う際に，

f

x及び

f

tを求める処理の終了を担保されている必要がある．並列計算できるデータを考慮した上で，図

3.3.7

^のように

stream

^を

0

^から

4

^の

5

つ用いて計算を行うようプログラムを書き換える．このときのオプティカルフロー計算のプログラムコードを図

3.3.8

^に示す．

cudaMemcpyAsync

関数で非同期で

CPU

^と

GPU

のメモリ間でデータ転送をする．

cudaStreamSynchronize

^{関数で指定} したストリームの処理が終わるのを待つ．この

2

つの関数を用いることで，データ依存性を壊さず並列で計算を進める．

stream 0 stream 1 stream 2 stream 3 stream 4

time

f x

f y

f t

S xx

S yy

S xy

S xt

S yt

S xx S yy

S xy S xy

S xt S yy

S yt S xx

S xt S xy

S yt S xy

u J DET

v

図

3.3.7:

マルチストリーム実装時のタイムライン

3.3.5 CPU

のハイブリッド並列の利用

§3.3.4

^{の実装では，}

GPU1

つ搭載のマシンでしか実行ができない．先行研究

[8]

^では

16

^台の計算機を相互結合し，

OpenMP

^と

MPI

のハイブリッド並列にて高速化を実現する．そこで，同じく

MPI

^{でノード間並列，}

OpenMP

でスレッド間並列を実装する．

§ 3.3.4

^{の実装では，}

CPU1

^コア分が

GPU

の制御を行っている．複数

CPU

1

単位として，早い者順でタスクを割り当てている．

1 cudaMemcpyAsync ( img1 , h_img1 , IMAGE , cudaMemcpyHostToDevice , stream [0]); // 非同期データ転送 2 cudaMemcpyAsync ( img2 , h_img2 , IMAGE , cudaMemcpyHostToDevice , stream [2]);

4 cudaStreamSynchronize ( stream [0]); // ストリーム同期 (データ転送の終了を待つ) 5

6 // stream 0 : m1 * dx -> fx

7 conv_kernel_shared <<< block , thread , shared1 , stream [0] > > >( img1 , fx , dx , w , h , f1 );

8 // stream 1 : m1 * dy -> fy

9 conv_kernel_shared <<< block , thread , shared1 , stream [1] > > >( img1 , fy , dy , w , h , f1 );

10 // stream 2 : m2 * gs -> ft

11 conv_kernel_shared <<< block , thread , shared1 , stream [2] > > >( img2 , ft , gs , w , h , f1 );

13 // stream 0 : fx * rf -> S_xx -> S_xx + ep

14 cudaStreamSynchronize ( stream [0]); // fx の計算が完了

15 mul_matrix_cuda <<< block , thread ,0 , stream [0] > > >( fx , fx , fxfx , w , h );

16 conv_kernel_shared <<< block , thread , shared2 , stream [0] > > >( fxfx , S_xx , rf , w , h , f2 );

17 add_matrix_cuda <<< block , thread ,0 , stream [0] > > >( S_xx , epsilon , w , h );

19 // stream 1 : fy * rf -> S_yy -> S_yy + ep

20 cudaStreamSynchronize ( stream [1]); // fy の計算が完了

21 mul_matrix_cuda <<< block , thread ,0 , stream [1] > > >( fy , fy , fyfy , w , h );

22 conv_kernel_shared <<< block , thread , shared2 , stream [1] > > >( fyfy , S_yy , rf , w , h , f2 );

23 add_matrix_cuda <<< block , thread ,0 , stream [1] > > >( S_yy , epsilon , w , h );

25 // stream 2 : fx , fy * rf -> S_xy

26 cudaStreamSynchronize ( stream [2]); // ft の計算が完了

27 mul_matrix_cuda <<< block , thread , 0, stream [2] > > >( fx , fy , fxfy , w , h );

28 conv_kernel_shared <<< block , thread , shared2 , stream [2] > > >( fxfy , S_xy , rf , w , h , f2 );

30 // stream 3 : fx , ft * rf -> S_xt

31 mul_matrix_cuda <<< block , thread , 0, stream [3] > > >( fx , ft , fxft , w , h );

32 conv_kernel_shared <<< block , thread , shared2 , stream [3] > > >( fxft , S_xt , rf , w , h , f2 );

34 // stream 4 : fy , ft * rf -> S_yt

35 mul_matrix_cuda <<< block , thread , 0, stream [4] > > >( fy , ft , fyft , w , h );

36 conv_kernel_shared <<< block , thread , shared2 , stream [4] > > >( fyft , S_yt , rf , w , h , f2 );

38 // stream 0 ,1

39 cudaStreamSynchronize ( stream [0]); // S_xx の計算が完了 40 cudaStreamSynchronize ( stream [1]); // S_yy の計算が完了

41 mul_matrix_cuda <<< block , thread ,0 , stream [1] > > >( S_xx , S_yy , S_xxyy , w , h ); // S_xx x S_yy 42 // stream 2

43 cudaStreamSynchronize ( stream [2]); // S_xy の計算が完了

44 mul_matrix_cuda <<< block , thread ,0 , stream [2] > > >( S_xy , S_xy , S_xyxy , w , h ); // S_xy x S_xy 45 // stream 3

46 cudaStreamSynchronize ( stream [3]); // S_xt の計算が完了

47 mul_matrix_cuda <<< block , thread ,0 , stream [3] > > >( S_xt , S_yy , S_xtyy , w , h ); // S_xt x S_yy 48 mul_matrix_cuda <<< block , thread ,0 , stream [3] > > >( S_xt , S_xy , S_xtxy , w , h ); // S_xt x S_xy 49 // stream 4

50 cudaStreamSynchronize ( stream [4]); // S_yt の計算が完了

51 mul_matrix_cuda <<< block , thread , 0, stream [4] > > >( S_yt , S_xx , S_ytxx , w , h ); // S_yt x S_xx 52 mul_matrix_cuda <<< block , thread , 0, stream [4] > > >( S_yt , S_xy , S_ytxy , w , h ); // S_yt x S_xy 53

54 // stream 1 ,2

55 cudaStreamSynchronize ( stream [1]); // S_xxyy の計算が完了 56 cudaStreamSynchronize ( stream [2]); // S_xyxy の計算が完了

57 sub_matrix_cuda <<< block , thread , 0, stream [0] > > >( S_xxyy , S_xyxy , J_DET , w , h );

ノード1

処理パターン数を計算 ← 全パターン数 / 全プロセス数

初期化

各画素パターンについてループ

結果の書き出し オプティカルフロー推定

錯視度合いの算出 円環状画像生成

スレッド1 3 4

オプティカルフロー推定 錯視度合いの算出

円環状画像生成 スレッド2

3 4 プロセスの起動

計算結果を親ノードに集約 プロセスの終了

MPIでプロセス並列

OpenMPでスレッド並列

GPU CPU

時間

図

3.3.9:

シミュレーションの流れ

(MPI+OpenMP+GPGPU)

1 # include <omp .h >

2 ^{# include} <mpi .h >

4 int world_size , my_rank ; 5 MPI_Init (& argc , & argv );

6 MPI_Comm_rank ( MPI_COMM_WORLD , & my_rank );

7 MPI_Comm_size ( MPI_COMM_WORLD , & world_size );

9 float res [ AllPattern ];

11 # pragma omp parallel for private ( PatternNumber ) schedule ( dynamic ,4) 12 for( PatternNumber = 0; PatternNumber < AllPattern ; PatternNumber ++){

13 if( omp_get_thread_num () == 0){ // スレッド番号0ならばG P Uで計算

14 ---- G P Uでの計算内容

----15 res [ PatternNumber ] = result_calc ; 16 } else { // スレッド番号1以上ならばC P Uで計算

17 C P Uでの計算内容

18 res [ PatternNumber ] = result_calc ; 19 }

21 // 結果の集約

22 if( my_rank != 0){ // 子ノードならば

23 ---- M P I _ S e n d関数で親ノードに結果を送信 ----24 } else {

25 for(int i =1; i < world_size ; i ++) {

26 ---- M P I _ S e n d関数で子ノードから結果を受信

----27 }

28 }

図

3.3.10:

シミュレーションのプログラムコードの一部

(MPI + OpenMP + GPGPU)

3.4 ^性能評価

3.4.1

各実装のシミュレーション時間の結果と評価

CPU (Intel Core i5-6500)

と

GPU (NVIDIA GeForce GTX 660)

を搭載したマシンにて，

§ 3.3.1

〜

3.3.4

の

GPGPU

における視覚数理モデルシミュレーションの高速化実装の計算時間を測定した．より詳

細なハードウェア・ソフトウェア環境は表

3.1

の通りである．なお，本来であれば

8

⁸

= 16 , 777 , 216

通りの入力パターンに対して計算を行うが，ここでは実験時間短縮のため，

256

^分の

1

^{通りにあた} る

4

⁸

= 65, 536

GPU

での計算に時間を要しているのはもちろんだが，

全体の約

20.6%

を

GPU

でのメモリ管理に，約

21.5%

を

CPU

と

GPU

(C)

メモリ確保及び解放の一括化をした実装では，メモリ管理時間をほぼ

0

^{秒まで削減できた．}

しかし，

GPU

での処理時間が増えた．これは

GPU

のカーネル内でポインタの付け替え等の処理が増えたことに起因していると考えられる．

CPU

^で

OpenMP

^を用いた

CPU

^{の計算時間と比べる}

と約

40.5%

削減できており，

GPGPU

による高い並列性を活かすことができた．

(D)

マルチストリーム実装では，メモリ管理を削減した実装と比べて計算時間は約

3.2%

しか削減できなかった．この実装時の

NVIDIA Visual Profiler

のタイムライン表示を図

3.4.2

^に示す．

CPU

と

GPU

間のデータ転送と畳み込み演算はオーバーラップできているが，畳み込み演算自体は並列に実行できていない．畳み込み演算に対して

GPU

の性能が足りていないことから，十分に並列実行できなかったと推測できる．

表

3.1:

ハードウェア・ソフトウェア環境

(GPU

^搭載，

1

^ノード

)

プロセッサ

Intel Core i5-6500

周波数

3.20 GHz

コア数

4

スレッド数

4

キャッシュ

6 MB

メモリ

32 GB

GPU NVIDIA GeForce GTX 660

CUDA

^コア数

960

プロセッサコアの周波数

980 MHz

メモリ

2048 MB

メモリ帯域幅

144.2 GB / sec

システム

IF PCI Express x16 Generation 3

OS Linux Mint 18.3

CXX compiler GNU 5.4.0

MPI

ライブラリ

OpenMPI 2.0.1

コンピュータビジョンライブラリ

OpenCV 3.2.0

CUDA C CUDA 9.0

20m23s

14m15s

8m53s 8m36s

(A)Prototype (B) Reducing

data transfer (C)

Centralized

memory management

(D)Multistream 0

500 1000 1500

Simulation Time [sec]

Processing(GPU)

Memory Management(GPU) Data transfer(CPU-GPU) Processing(CPU) Sum(CPU+GPU)

図

3.4.1: 65,536

通りのシミュレーション時間と割合

(

単一ノード

)

OpticalFlow Calculation for each image

Data transfer

CPU→GPU Data transfer

GPU→CPU

Convolution

Multi-stream

Overlap

time time

OpticalFlow Calculation for each image

Data transfer

CPU→GPU Data transfer

GPU→CPU

Convolution

Single-stream

Before Multi-streaming

After Multi-streaming

図

3.4.2:

マルチストリーム化実装前後の

NVIDIA Visual Profiler

^{のタイムライン比較}

3.4.2 GPU

クラスタを用いたシミュレーション時間の結果と評価

CPU (Intel Core i5-8400)

^と

GPU (NVIDIA GeForce GTX 660)

^{を搭載した}

16

^台の

GPU

^クラスタ

にて，

§ 3.3.5

^の

GPGPU

における視覚数理モデルシミュレーションの高速化実装の計算時間を測定

した．より詳細なハードウェア・ソフトウェア環境は表

3.2

CPU 6

コアで並列処理されている．

(E)

^{と比べると}

CPU

GPU

1

台は

CPU 4.5

コア分の計算能力であると予測したことから，

CPU

及び

GPU

の全リソースを利用

する実装では

CPU 5 + 4.5 = 9.5

コア分の計算能力であると推測できる．実際に

CPU 1

^{コアのみで} ある

(E)

^{と比べて約}

8.7

^{倍速であり，}

OpenMP

を用いたことによるオーバーヘッドはあるものの十分に計算機リソースが利用できているといえる．

また，先行研究

[3]

^{の視覚数理モデルは}

MATLAB

にて実装されており，元々は

6

^{日ほどかかる} 計算を本研究では

CUDA

^で記述し

GPU

クラスタを用いることで

2%

^{以下に削減できた．}

ドキュメント内視覚数理モデルシミュレーションの並列化と錯視画像の探索 (ページ 76-115)