Mereo

ループ展開でC/C++畳み込み処理の高速化

目次

  1. 1. はじめに
  2. 2. 概要説明
  3. 3. 畳み込み処理
  4. 4. 今回の畳み込み仕様
  5. 5. 動作環境など、条件
  6. 6. 初期実装
  7. 7. 高速化
    1. 7.1. ループまわりを少し
    2. 7.2. さらにループ展開
    3. 7.3. おわりに

はじめに

呉高専 Advent Calendar 2018
18日目の記事です。

本来は「Python vs C++ 高速化対決」というタイトルでしたが、
諸事情でボツになりました。

概要説明

この記事ではC/C++で畳み込み処理を実装し、ループ展開を中心に高速化します。

(今回のコードは全体的にかなり適当です。言語哲学者の方はお帰り下さい)

畳み込み処理

みんな知ってる画像処理の分野では必須な処理です。

どんな処理なのかを知りたい方は
ここ
などを参考にしてください。

今回の畳み込み仕様

  • 画像サイズ：720x480
  • フィルタサイズ：5x5
  • padding：2　(右と下に4)

動作環境など、条件

  • CPU：インテル® Core™ i7-4770 プロセッサー
    • ベース：3.40 GHz
    • コア：4
    • キャッシュ：8 MB
    • 拡張命令：AVX2対応
  • コンパイラ：VC++ v141
  • 最適化オプション /O2

初期実装

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
void Convolution2d(UCHAR * src, UCHAR * dst, UCHAR * k, UINT height, UINT width) {
	
	unsigned short sum = 0;
	for (UINT i = 0; i < height; i++) {
		for (UINT j = 0; j < width; j++) {
			for (UINT k_i = 0; k_i < 5; k_i++) {
				for (UINT k_j = 0; k_j < 5; k_j++) {
					sum += (unsigned short)src[(i*(width+4) + j) + (k_i * (width+4) + k_j)] * (unsigned short)k[k_i*5+k_j];
				}
			}
			dst[i*(width) + j] = (UCHAR)(sum/25);
			sum = 0;
		}
	}
}
1
Time: 5.09439[msec]

愚直に組んだ4重ループでの処理ですね。頭悪そうです。

一応、画像データとフィルタはあらかじめ32Byteでアラインしてあります。
(SIMDは使う予定ないですが。)

計測時間には10000回で計測し1回分の時間を用いています。

さあ、高速化していきましょう。

高速化

ループまわりを少し

やることと理由

  • 二つほど完全にループ展開(邪魔なので、条件文通過回数が減るので)
  • アドレスを使用したループ制御(配列添え字のインデックス計算が無駄なので、インクリメントが無駄なので)
  • フィルタ2つ分を同ループ内で処理(少しループ展開になる、メモリ最適化にもなる)
  • forからwhileに変更(forは遅いので)

変更済みコードは以下の通りです。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
void Convolution2d_opt(UCHAR * src, UCHAR * dst, UCHAR * k, UINT height, UINT width) {
	UCHAR _kx[50];
	UCHAR* kx = _kx;
	for (int i = 0; i < 25; i++) {
		kx[i * 2] = k[i];
		kx[i * 2+1] = k[i];
	}

	UCHAR* kx_start = kx;
	UCHAR* src_start = src;


	USHORT sum0 = 0, sum1 = 0;

	UINT width4 = width + 4, height4 = height + 4;
	UCHAR* src_i_end = src + height * width4;
	UCHAR* src_j_end;
	UCHAR * src_i = src_start , *src_j;
	UCHAR* src_cur;

	while (src_i < src_i_end) {
		src_j_end = src_i + width;
		src_j = src_i;
		while (src_j < src_j_end) {
			src_cur = src_j;
			kx = kx_start;

			sum0 += (USHORT)*(src_cur + 0) * (USHORT)*(kx + 0);
			sum1 += (USHORT)*(src_cur + 1) * (USHORT)*(kx + 1);
			sum0 += (USHORT)*(src_cur + 1) * (USHORT)*(kx + 2);
			sum1 += (USHORT)*(src_cur + 2) * (USHORT)*(kx + 3);
			sum0 += (USHORT)*(src_cur + 2) * (USHORT)*(kx + 4);
			sum1 += (USHORT)*(src_cur + 3) * (USHORT)*(kx + 5);
			sum0 += (USHORT)*(src_cur + 3) * (USHORT)*(kx + 6);
			sum1 += (USHORT)*(src_cur + 4) * (USHORT)*(kx + 7);
			sum0 += (USHORT)*(src_cur + 4) * (USHORT)*(kx + 8);
			sum1 += (USHORT)*(src_cur + 5) * (USHORT)*(kx + 9);

			kx += 10; src_cur += width4;

			sum0 += (USHORT)*(src_cur + 0) * (USHORT)*(kx + 0);
			sum1 += (USHORT)*(src_cur + 1) * (USHORT)*(kx + 1);
			sum0 += (USHORT)*(src_cur + 1) * (USHORT)*(kx + 2);
			sum1 += (USHORT)*(src_cur + 2) * (USHORT)*(kx + 3);
			sum0 += (USHORT)*(src_cur + 2) * (USHORT)*(kx + 4);
			sum1 += (USHORT)*(src_cur + 3) * (USHORT)*(kx + 5);
			sum0 += (USHORT)*(src_cur + 3) * (USHORT)*(kx + 6);
			sum1 += (USHORT)*(src_cur + 4) * (USHORT)*(kx + 7);
			sum0 += (USHORT)*(src_cur + 4) * (USHORT)*(kx + 8);
			sum1 += (USHORT)*(src_cur + 5) * (USHORT)*(kx + 9);

			kx += 10; src_cur += width4;

			sum0 += (USHORT)*(src_cur + 0) * (USHORT)*(kx + 0);
			sum1 += (USHORT)*(src_cur + 1) * (USHORT)*(kx + 1);
			sum0 += (USHORT)*(src_cur + 1) * (USHORT)*(kx + 2);
			sum1 += (USHORT)*(src_cur + 2) * (USHORT)*(kx + 3);
			sum0 += (USHORT)*(src_cur + 2) * (USHORT)*(kx + 4);
			sum1 += (USHORT)*(src_cur + 3) * (USHORT)*(kx + 5);
			sum0 += (USHORT)*(src_cur + 3) * (USHORT)*(kx + 6);
			sum1 += (USHORT)*(src_cur + 4) * (USHORT)*(kx + 7);
			sum0 += (USHORT)*(src_cur + 4) * (USHORT)*(kx + 8);
			sum1 += (USHORT)*(src_cur + 5) * (USHORT)*(kx + 9);

			kx += 10; src_cur += width4;

			sum0 += (USHORT)*(src_cur + 0) * (USHORT)*(kx + 0);
			sum1 += (USHORT)*(src_cur + 1) * (USHORT)*(kx + 1);
			sum0 += (USHORT)*(src_cur + 1) * (USHORT)*(kx + 2);
			sum1 += (USHORT)*(src_cur + 2) * (USHORT)*(kx + 3);
			sum0 += (USHORT)*(src_cur + 2) * (USHORT)*(kx + 4);
			sum1 += (USHORT)*(src_cur + 3) * (USHORT)*(kx + 5);
			sum0 += (USHORT)*(src_cur + 3) * (USHORT)*(kx + 6);
			sum1 += (USHORT)*(src_cur + 4) * (USHORT)*(kx + 7);
			sum0 += (USHORT)*(src_cur + 4) * (USHORT)*(kx + 8);
			sum1 += (USHORT)*(src_cur + 5) * (USHORT)*(kx + 9);


			kx += 10; src_cur += width4;

			sum0 += (USHORT)*(src_cur + 0) * (USHORT)*(kx + 0);
			sum1 += (USHORT)*(src_cur + 1) * (USHORT)*(kx + 1);
			sum0 += (USHORT)*(src_cur + 1) * (USHORT)*(kx + 2);
			sum1 += (USHORT)*(src_cur + 2) * (USHORT)*(kx + 3);
			sum0 += (USHORT)*(src_cur + 2) * (USHORT)*(kx + 4);
			sum1 += (USHORT)*(src_cur + 3) * (USHORT)*(kx + 5);
			sum0 += (USHORT)*(src_cur + 3) * (USHORT)*(kx + 6);
			sum1 += (USHORT)*(src_cur + 4) * (USHORT)*(kx + 7);
			sum0 += (USHORT)*(src_cur + 4) * (USHORT)*(kx + 8);
			sum1 += (USHORT)*(src_cur + 5) * (USHORT)*(kx + 9);


			*dst = sum0 / 25; dst++;
			*dst = sum1 / 25; dst++;
			sum0 = 0;
			sum1 = 0;

			src_j += 2;
		}
		src_i += width4;
	}
}
1
Time: 2.98714[msec]

約1.71倍になりました。

いつもならここで妥協しますが、ここからさらにループ展開でごり押します。

さらにループ展開

このままループ展開したら可読性が低下します。
defineを使いましょう。

やることと理由

  • フィルタ4つ分を同ループ内で処理(少しループ展開になる、メモリ最適化にもなる)
  • define使ってごり押しループ展開

変更済みコードは以下の通りです。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
#define Calc_1(level)\
*(sum + 0) += (UCHAR)*(src_cur + 0 + level) * (UCHAR)*(kx + level*4 + 0);\
*(sum + 1) += (UCHAR)*(src_cur + 1 + level) * (UCHAR)*(kx + level*4 + 1);\
*(sum + 2) += (UCHAR)*(src_cur + 2 + level) * (UCHAR)*(kx + level*4 + 2);\
*(sum + 3) += (UCHAR)*(src_cur + 3 + level) * (UCHAR)*(kx + level*4 + 3)

#define Calc_5()\
Calc_1(0);\
Calc_1(1);\
Calc_1(2);\
Calc_1(3);\
Calc_1(4)

#define Calc_25()\
src_cur = src_j;\
kx = kx_start;\
Calc_5();\
kx += 20; src_cur += width4;\
Calc_5();\
kx += 20; src_cur += width4;\
Calc_5();\
kx += 20; src_cur += width4;\
Calc_5();\
kx += 20; src_cur += width4;\
Calc_5();\
*dst = *(sum + 0) / 25; dst++;\
*dst = *(sum + 1) / 25; dst++;\
*dst = *(sum + 2) / 25; dst++;\
*dst = *(sum + 3) / 25; dst++;\
*(sum + 0) = 0;\
*(sum + 1) = 0;\
*(sum + 2) = 0;\
*(sum + 3) = 0;



void Convolution2d_opt(UCHAR * src, UCHAR * dst, UCHAR * k, UINT height, UINT width) {
	UCHAR _kx[100];
	UCHAR* kx = _kx;
	UCHAR* kx_start = kx;
	for (int i = 0; i < 25; i++) {
		kx[i * 4] = k[i];
		kx[i * 4 + 1] = k[i];
		kx[i * 4 + 2] = k[i];
		kx[i * 4 + 3] = k[i];
	}

	UCHAR* src_start = src;


	USHORT sum[4] = { 0,0,0,0 };

	UINT width4 = width + 4, height4 = height + 4;
	UCHAR* src_i_end = src + height * width4;
	UCHAR* src_j_end;
	UCHAR * src_i = src_start , *src_j;
	UCHAR* src_cur;

	while (src_i < src_i_end) {
		src_j_end = src_i + width;
		src_j = src_i;
		while (src_j < src_j_end) {
			Calc_25();
			src_j += 4;
			Calc_25();
			src_j += 4;
			Calc_25();
			src_j += 4;
			Calc_25();
			src_j += 4;
			Calc_25();
			src_j += 4;
			Calc_25();
			src_j += 4;
			Calc_25();
			src_j += 4;
			Calc_25();
			src_j += 4;
			Calc_25();
			src_j += 4;
			Calc_25();
			src_j += 4;
		}
		src_i += width4;
	}
}
1
Time: 2.73656[msec]

約1.86倍まで上がりました。
ループ展開は素晴らしいですね！

しかし、これ以上ループ展開すると性能が悪化してしまいます。
コードが長すぎてコンパイラが最適化をあきらめたのでしょう。

何事にも限度というものがありますね。ここまでにしておきましょう。

おわりに

今回はループまわりに絞ってちょっと高速化してみました。

ごり押しループ展開で最後の悪あがき程度にはなることが分かってもらえたと思います。

しかしループ展開にも「コード量の肥大化」「コンパイル時間が増える」等のデメリットも存在します。

用法用量にはくれぐれもご注意を。

因みに、他にも最適化できる場所はまだまあります。

最速を目指したい方は、「SIMD、キャッシュ、メモリ最適化、GPGPU」
等の単語で調べてみるといいでしょう。