Post

Convolution with SIMD

Posted Updated
By
20 min read
Convolution with SIMD

Convolution with SIMD

Prerequisites

1

SIMD

1. Naive Code

About 5119 by 5119 image.

Processing time: 111ms

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

uint8_t* pU8SrcBuf = (uint8_t*)imgSrc.data;
uint8_t* pU8DstBuf = (uint8_t*)imgDst.data;

for(int i32Row = 2; i32Row < imgSrc.rows - 2; ++i32Row)
{
    uint8_t* pU8SrcBufRow = pU8SrcBuf + i32Row * imgSrc.step + 2;
    uint8_t* pU8DstBufRow = pU8DstBuf + i32Row * imgDst.step + 2;

    for(int i32Col = 2; i32Col < imgSrc.cols - 2; ++i32Col)
    {
        int i32Sum = 0;

        for(int i = -2; i <= 2; ++i)
        {
            for(int j = -2; j <= 2; ++j)
                i32Sum += *(pU8SrcBufRow + j + (i * imgSrc.step));
        }

        *pU8DstBufRow = uint8_t(i32Sum / 25);

        ++pU8SrcBufRow;
        ++pU8DstBufRow;
    }
}

Optimization Points

Cache locality

approach Contiguous memory about row, not column

1
2
3
4
5

for(int i = -2; i <= 2; ++i)
{
    for(int j = -2; j <= 2; ++j)
        i32Sum += *(pU8SrcBufRow + j + (i * imgSrc.step));
}
Don’t make double point dereference cost
1
2

uint8_t* pU8SrcBufRow = pU8SrcBuf + i32Row * imgSrc.step + 2;
uint8_t* pU8DstBufRow = pU8DstBuf + i32Row * imgDst.step + 2;
Don’t make duplicate cost
1

*pU8DstBufRow = uint8_t(i32Sum / 25);

or

1

*pU8DstBufRow = (i32Sum * 2621) >> 16;

Because division cost is very higher than mul and shift. But it is approximate value.

2. Naive Code Advanced

Processing time: 85ms

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

uint8_t* pU8SrcBuf = (uint8_t*)imgSrc.data;
uint8_t* pU8DstBuf = (uint8_t*)imgDst.data;

for(int i32Row = 2; i32Row < imgSrc.rows - 2; ++i32Row)
{
    uint8_t* pU8SrcBufRow0 = pU8SrcBuf + (i32Row - 2) * imgSrc.step + 2;
    uint8_t* pU8SrcBufRow1 = pU8SrcBuf + (i32Row - 1) * imgSrc.step + 2;
    uint8_t* pU8SrcBufRow2 = pU8SrcBuf + i32Row * imgSrc.step + 2;
    uint8_t* pU8SrcBufRow3 = pU8SrcBuf + (i32Row + 1) * imgSrc.step + 2;
    uint8_t* pU8SrcBufRow4 = pU8SrcBuf + (i32Row + 2) * imgSrc.step + 2;

    uint8_t* pU8DstBufRow = pU8DstBuf + i32Row * imgDst.step;

    for(int i32Col = 2; i32Col < imgSrc.cols - 2; ++i32Col)
    {
        int i32Sum = 0;

        i32Sum += *(pU8SrcBufRow0 - 2) + *(pU8SrcBufRow0 - 1) + *(pU8SrcBufRow0) + *(pU8SrcBufRow0 + 1) + *(pU8SrcBufRow0 + 2);
        i32Sum += *(pU8SrcBufRow1 - 2) + *(pU8SrcBufRow1 - 1) + *(pU8SrcBufRow1) + *(pU8SrcBufRow1 + 1) + *(pU8SrcBufRow1 + 2);
        i32Sum += *(pU8SrcBufRow2 - 2) + *(pU8SrcBufRow2 - 1) + *(pU8SrcBufRow2) + *(pU8SrcBufRow2 + 1) + *(pU8SrcBufRow2 + 2);
        i32Sum += *(pU8SrcBufRow3 - 2) + *(pU8SrcBufRow3 - 1) + *(pU8SrcBufRow3) + *(pU8SrcBufRow3 + 1) + *(pU8SrcBufRow3 + 2);
        i32Sum += *(pU8SrcBufRow4 - 2) + *(pU8SrcBufRow4 - 1) + *(pU8SrcBufRow4) + *(pU8SrcBufRow4 + 1) + *(pU8SrcBufRow4 + 2);

        *pU8DstBufRow = uint8_t(i32Sum / 25);

        pU8SrcBufRow0++;
        pU8SrcBufRow1++;
        pU8SrcBufRow2++;
        pU8SrcBufRow3++;
        pU8SrcBufRow4++;
        pU8DstBufRow++;
    }
}

Optimization Points

Loop unrolling
1
2
3
4
5

i32Sum += *(pU8SrcBufRow0 - 2) + *(pU8SrcBufRow0 - 1) + *(pU8SrcBufRow0) + *(pU8SrcBufRow0 + 1) + *(pU8SrcBufRow0 + 2);
i32Sum += *(pU8SrcBufRow1 - 2) + *(pU8SrcBufRow1 - 1) + *(pU8SrcBufRow1) + *(pU8SrcBufRow1 + 1) + *(pU8SrcBufRow1 + 2);
i32Sum += *(pU8SrcBufRow2 - 2) + *(pU8SrcBufRow2 - 1) + *(pU8SrcBufRow2) + *(pU8SrcBufRow2 + 1) + *(pU8SrcBufRow2 + 2);
i32Sum += *(pU8SrcBufRow3 - 2) + *(pU8SrcBufRow3 - 1) + *(pU8SrcBufRow3) + *(pU8SrcBufRow3 + 1) + *(pU8SrcBufRow3 + 2);
i32Sum += *(pU8SrcBufRow4 - 2) + *(pU8SrcBufRow4 - 1) + *(pU8SrcBufRow4) + *(pU8SrcBufRow4 + 1) + *(pU8SrcBufRow4 + 2);

3. SIMD

Processing time: 12ms

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

uint8_t* pU8SrcBuf = (uint8_t*)imgSrc.data;
uint8_t* pU8DstBuf = (uint8_t*)imgDst.data;

__m256i m256iKernelDiv = _mm256_set1_epi16(25);
__m256i m256iZero = _mm256_setzero_si256();

    for(int i32Row = 2; i32Row < imgSrc.rows - 2; ++i32Row)
{
    uint8_t* pU8SrcBufRow = pU8SrcBuf + i32Row * imgSrc.step + 2;
    uint8_t* pU8DstBufRow = pU8DstBuf + i32Row * imgDst.step + 2;

    int i32Col = 2;

    for(; i32Col + 32 < imgSrc.cols - 2; i32Col +=32)
    {
        __m256i m256iSumlo = _mm256_setzero_si256();
        __m256i m256iSumhi = _mm256_setzero_si256();

        for(int i = -2; i <= 2; ++i)
        {
            for(int j = -2; j <= 2; ++j)
            {
                __m256i m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + j + (i * imgSrc.step)));
                __m256i m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero); // latency 1
                __m256i m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero); // latency 1

                m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
                m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);
            }
        }

        m256iSumlo = _mm256_div_epi16(m256iSumlo, m256iKernelDiv); // latency, maybe very high
        m256iSumhi = _mm256_div_epi16(m256iSumhi, m256iKernelDiv);
        __m256i m256iDst = _mm256_packus_epi16(m256iSumlo, m256iSumhi); // latency 3

        _mm256_storeu_si256((__m256i*)pU8DstBufRow, m256iDst);
        pU8SrcBufRow += 32; 
        pU8DstBufRow += 32; 
    }

    for(; i32Col < imgSrc.cols - 2; ++i32Col)
    {
        int i32Sum = 0;

        for(int i = -2; i <= 2; ++i)
        {
            for(int j = -2; j <= 2; ++j)
                i32Sum += *(pU8SrcBufRow + j + (i * imgSrc.step));
        }

        *pU8DstBufRow = uint8_t(i32Sum / 25);

        ++pU8SrcBufRow;
        ++pU8DstBufRow;
    }
}

Optimization Points

SIMD

Convolution is to be sumed from multiple data. considering overflow and seperate the memory, process and pack again.

1
2
3
4
5
6
7
8

_mm256_setzero_si256()  // latency 1
_mm256_loadu_si256()    // latency 7
_mm256_unpacklo_epi8()  // latency 1
_mm256_unpackhi_epi8()  // latency 1
_mm256_add_epi16()      // latency 1
_mm256_div_epi16()      // ?
_mm256_packus_epi16()   // latency 3
_mm256_storeu_si256()   // latency 1
  • _mm256_setzero_si256
1

dst[MAX:0] := 0
  • _mm256_loadu_si256()
1
2

dst[255:0] := MEM[mem_addr+255:mem_addr]
dst[MAX:256] := 0
  • _mm256_unpacklo_epi8()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

DEFINE INTERLEAVE_BYTES(src1[127:0], src2[127:0]) {
	dst[7:0] := src1[7:0] 
	dst[15:8] := src2[7:0] 
	dst[23:16] := src1[15:8] 
	dst[31:24] := src2[15:8] 
	dst[39:32] := src1[23:16] 
	dst[47:40] := src2[23:16] 
	dst[55:48] := src1[31:24] 
	dst[63:56] := src2[31:24] 
	dst[71:64] := src1[39:32]
	dst[79:72] := src2[39:32] 
	dst[87:80] := src1[47:40] 
	dst[95:88] := src2[47:40] 
	dst[103:96] := src1[55:48] 
	dst[111:104] := src2[55:48] 
	dst[119:112] := src1[63:56] 
	dst[127:120] := src2[63:56] 
	RETURN dst[127:0]
}
dst[127:0] := INTERLEAVE_BYTES(a[127:0], b[127:0])
dst[255:128] := INTERLEAVE_BYTES(a[255:128], b[255:128])
dst[MAX:256] := 0
  • _mm256_unpackhi_epi8()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

DEFINE INTERLEAVE_HIGH_BYTES(src1[127:0], src2[127:0]) {
	dst[7:0] := src1[71:64] 
	dst[15:8] := src2[71:64] 
	dst[23:16] := src1[79:72] 
	dst[31:24] := src2[79:72] 
	dst[39:32] := src1[87:80] 
	dst[47:40] := src2[87:80] 
	dst[55:48] := src1[95:88] 
	dst[63:56] := src2[95:88] 
	dst[71:64] := src1[103:96] 
	dst[79:72] := src2[103:96] 
	dst[87:80] := src1[111:104] 
	dst[95:88] := src2[111:104] 
	dst[103:96] := src1[119:112] 
	dst[111:104] := src2[119:112] 
	dst[119:112] := src1[127:120] 
	dst[127:120] := src2[127:120] 
	RETURN dst[127:0]	
}
dst[127:0] := INTERLEAVE_HIGH_BYTES(a[127:0], b[127:0])
dst[255:128] := INTERLEAVE_HIGH_BYTES(a[255:128], b[255:128])
dst[MAX:256] := 0
  • _mm256_add_epi16()
1
2
3
4
5

FOR j := 0 to 15
	i := j*16
	dst[i+15:i] := a[i+15:i] + b[i+15:i]
ENDFOR
dst[MAX:256] := 0
  • _mm256_div_epi16()
1
2
3
4
5
6
7
8

FOR j := 0 to 15
	i := 16*j
	IF b[i+15:i] == 0
		#DE
	FI
	dst[i+15:i] := Truncate16(a[i+15:i] / b[i+15:i])
ENDFOR
dst[MAX:256] := 0
  • _mm256_packus_epi16()
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

dst[7:0] := SaturateU8(a[15:0])
dst[15:8] := SaturateU8(a[31:16])
dst[23:16] := SaturateU8(a[47:32])
dst[31:24] := SaturateU8(a[63:48])
dst[39:32] := SaturateU8(a[79:64])
dst[47:40] := SaturateU8(a[95:80])
dst[55:48] := SaturateU8(a[111:96])
dst[63:56] := SaturateU8(a[127:112])
dst[71:64] := SaturateU8(b[15:0])
dst[79:72] := SaturateU8(b[31:16])
dst[87:80] := SaturateU8(b[47:32])
dst[95:88] := SaturateU8(b[63:48])
dst[103:96] := SaturateU8(b[79:64])
dst[111:104] := SaturateU8(b[95:80])
dst[119:112] := SaturateU8(b[111:96])
dst[127:120] := SaturateU8(b[127:112])
dst[135:128] := SaturateU8(a[143:128])
dst[143:136] := SaturateU8(a[159:144])
dst[151:144] := SaturateU8(a[175:160])
dst[159:152] := SaturateU8(a[191:176])
dst[167:160] := SaturateU8(a[207:192])
dst[175:168] := SaturateU8(a[223:208])
dst[183:176] := SaturateU8(a[239:224])
dst[191:184] := SaturateU8(a[255:240])
dst[199:192] := SaturateU8(b[143:128])
dst[207:200] := SaturateU8(b[159:144])
dst[215:208] := SaturateU8(b[175:160])
dst[223:216] := SaturateU8(b[191:176])
dst[231:224] := SaturateU8(b[207:192])
dst[239:232] := SaturateU8(b[223:208])
dst[247:240] := SaturateU8(b[239:224])
dst[255:248] := SaturateU8(b[255:240])
dst[MAX:256] := 0
  • _mm256_storeu_si256()
1

MEM[mem_addr+255:mem_addr] := a[255:0]

4. SIMD with loop unrolling

Processing time: 16ms

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

__m256i m256iKernelDiv = _mm256_set1_epi16(25);
__m256i m256iZero = _mm256_setzero_si256();

for(int i32Row = 2; i32Row < imgSrc.rows - 2; ++i32Row)
{
    uint8_t* pU8SrcBufRow = pU8SrcBuf + i32Row * imgSrc.step + 2;
    uint8_t* pU8DstBufRow = pU8DstBuf + i32Row * imgDst.step + 2;

    int i32Col = 2;

    for(; i32Col + 32 < imgSrc.cols - 2; i32Col +=32)
    {
        __m256i m256iSumlo = _mm256_setzero_si256();
        __m256i m256iSumhi = _mm256_setzero_si256();
        __m256i m256iSrc, m256iSrclo, m256iSrchi; 

        for(int i = -2; i <= 2; ++i)
        {
            m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow - 2 + (i * imgSrc.step)));
            m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
            m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
            m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
            m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

            m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow - 1 + (i * imgSrc.step)));
            m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
            m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
            m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
            m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

            m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + (i * imgSrc.step)));
            m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
            m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
            m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
            m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

            m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + 1 + (i * imgSrc.step)));
            m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
            m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
            m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
            m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

            m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + 2 + (i * imgSrc.step)));
            m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
            m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
            m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
            m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);
        }

        m256iSumlo = _mm256_div_epi16(m256iSumlo, m256iKernelDiv);
        m256iSumhi = _mm256_div_epi16(m256iSumhi, m256iKernelDiv);
        __m256i m256iDst = _mm256_packus_epi16(m256iSumlo, m256iSumhi);

        _mm256_storeu_si256((__m256i*)pU8DstBufRow, m256iDst);
        pU8SrcBufRow += 32; 
        pU8DstBufRow += 32; 
    }

    for(; i32Col < imgSrc.cols - 2; ++i32Col)
    {
        int i32Sum = 0;

        for(int i = -2; i <= 2; ++i)
        {
            for(int j = -2; j <= 2; ++j)
                i32Sum += *(pU8SrcBufRow + j + (i * imgSrc.step));
        }

        *pU8DstBufRow = uint8_t(i32Sum / 25);

        ++pU8SrcBufRow;
        ++pU8DstBufRow;
    }
}

Optimization Points

Bad case

Loop unrolling is usually good but sometimes bad.

Just test!

5. SIMD with loop unrolling extremely

Processing time: 12ms

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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200

__m256i m256iKernelDiv = _mm256_set1_epi16(25);
__m256i m256iZero = _mm256_setzero_si256();

for(int i32Row = 2; i32Row < imgSrc.rows - 2; ++i32Row)
{
    uint8_t* pU8SrcBufRow = pU8SrcBuf + i32Row * imgSrc.step + 2;
    uint8_t* pU8DstBufRow = pU8DstBuf + i32Row * imgDst.step + 2;

    int i32Col = 2;

    for(; i32Col + 32 < imgSrc.cols - 2; i32Col +=32)
    {
        __m256i m256iSumlo = _mm256_setzero_si256();
        __m256i m256iSumhi = _mm256_setzero_si256();
        __m256i m256iSrc, m256iSrclo, m256iSrchi; 

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow - 2 + (-2 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow - 1 + (-2 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + (-2 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + 1 + (-2 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + 2 + (-2 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);


        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow - 2 + (-1 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow - 1 + (-1 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + (-1 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + 1 + (-1 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + 2 + (-1 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);



        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow - 2));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow - 1));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + 1));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + 2));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);



        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow - 2 + (1 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow - 1 + (1 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + (1 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + 1 + (1 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + 2 + (1 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);



        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow - 2 + (2 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow - 1 + (2 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + (2 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + 1 + (2 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);

        m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + 2 + (2 * imgSrc.step)));
        m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero);
        m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero);
        m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
        m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);



        m256iSumlo = _mm256_div_epi16(m256iSumlo, m256iKernelDiv);
        m256iSumhi = _mm256_div_epi16(m256iSumhi, m256iKernelDiv);
        __m256i m256iDst = _mm256_packus_epi16(m256iSumlo, m256iSumhi);

        _mm256_storeu_si256((__m256i*)pU8DstBufRow, m256iDst);
        pU8SrcBufRow += 32; 
        pU8DstBufRow += 32; 
    }

    for(; i32Col < imgSrc.cols - 2; ++i32Col)
    {
        int i32Sum = 0;

        for(int i = -2; i <= 2; ++i)
        {
            for(int j = -2; j <= 2; ++j)
                i32Sum += *(pU8SrcBufRow + j + (i * imgSrc.step));
        }

        *pU8DstBufRow = uint8_t(i32Sum / 25);

        ++pU8SrcBufRow;
        ++pU8DstBufRow;
    }
}

Optimization Points

similar case

Loop unrolling is usually good but sometimes bad.

Just test!

6. SIMD with replacing divide to shift

Processing time: 9.3ms

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

__m256i m256iKernelDivTrick = _mm256_set1_epi16(2621);
__m256i m256iZero = _mm256_setzero_si256();

for(int i32Row = 2; i32Row < imgSrc.rows - 2; ++i32Row)
{
    uint8_t* pU8SrcBufRow = pU8SrcBuf + i32Row * imgSrc.step + 2;
    uint8_t* pU8DstBufRow = pU8DstBuf + i32Row * imgDst.step + 2;

    int i32Col = 2;

    for(; i32Col + 32 < imgSrc.cols - 2; i32Col +=32)
    {
        __m256i m256iSumlo = _mm256_setzero_si256();
        __m256i m256iSumhi = _mm256_setzero_si256();

        for(int i = -2; i <= 2; ++i)
        {
            for(int j = -2; j <= 2; ++j)
            {
                __m256i m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRow + j + (i * imgSrc.step)));
                __m256i m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero); // latency 1
                __m256i m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero); // latency 1

                m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
                m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);
            }
        }

        // m256iSumlo = _mm256_div_epi16(m256iSumlo, m256iKernelDiv); // latency, maybe very high
        // m256iSumhi = _mm256_div_epi16(m256iSumhi, m256iKernelDiv);
        m256iSumlo = _mm256_mulhi_epu16(m256iSumlo, m256iKernelDivTrick); // latency 5
        m256iSumhi = _mm256_mulhi_epu16(m256iSumhi, m256iKernelDivTrick);

        __m256i m256iDst = _mm256_packus_epi16(m256iSumlo, m256iSumhi); // latency 3

        _mm256_storeu_si256((__m256i*)pU8DstBufRow, m256iDst);
        pU8SrcBufRow += 32; 
        pU8DstBufRow += 32; 
    }

    for(; i32Col < imgSrc.cols - 2; ++i32Col)
    {
        int i32Sum = 0;

        for(int i = -2; i <= 2; ++i)
        {
            for(int j = -2; j <= 2; ++j)
                i32Sum += *(pU8SrcBufRow + j + (i * imgSrc.step));
        }

        *pU8DstBufRow = uint8_t(i32Sum / 25);

        ++pU8SrcBufRow;
        ++pU8DstBufRow;
    }
}

Optimization Points

SIMD

Divide is high latency usually, So if we can we should replace division to other

1

_mm256_mulhi_epu16()  // latency 5
  • _mm256_mulhi_epu16()
1
2
3
4
5
6

FOR j := 0 to 15
	i := j*16
	tmp[31:0] := a[i+15:i] * b[i+15:i]
	dst[i+15:i] := tmp[31:16]
ENDFOR
dst[MAX:256] := 0

7. SIMD with reducing duplicate operation

Processing time: 8.5ms

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

__m256i m256iKernelDivTrick = _mm256_set1_epi16(2621);
__m256i m256iZero = _mm256_setzero_si256();

for(int i32Row = 2; i32Row < imgSrc.rows - 2; ++i32Row)
{
    uint8_t* pU8SrcBufRow = pU8SrcBuf + i32Row * imgSrc.step + 2;
    uint8_t* pU8DstBufRow = pU8DstBuf + i32Row * imgDst.step + 2;

    int i32Col = 2;

    for(; i32Col + 32 < imgSrc.cols - 2; i32Col +=32)
    {
        __m256i m256iSumlo = _mm256_setzero_si256();
        __m256i m256iSumhi = _mm256_setzero_si256();

        for(int i = -2; i <= 2; ++i)
        {
            uint8_t* pU8SrcBufRowCur = pU8SrcBufRow + i * imgSrc.step;

            for(int j = -2; j <= 2; ++j)
            {
                __m256i m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRowCur + j));
                __m256i m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero); // latency 1
                __m256i m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero); // latency 1

                m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
                m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);
            }
        }

        m256iSumlo = _mm256_mulhi_epu16(m256iSumlo, m256iKernelDivTrick); // latency 5
        m256iSumhi = _mm256_mulhi_epu16(m256iSumhi, m256iKernelDivTrick);

        __m256i m256iDst = _mm256_packus_epi16(m256iSumlo, m256iSumhi); // latency 3

        _mm256_storeu_si256((__m256i*)pU8DstBufRow, m256iDst);
        pU8SrcBufRow += 32; 
        pU8DstBufRow += 32; 
    }

    for(; i32Col < imgSrc.cols - 2; ++i32Col)
    {
        int i32Sum = 0;

        for(int i = -2; i <= 2; ++i)
        {
            for(int j = -2; j <= 2; ++j)
                i32Sum += *(pU8SrcBufRow + j + (i * imgSrc.step));
        }

        *pU8DstBufRow = uint8_t(i32Sum / 25);

        ++pU8SrcBufRow;
        ++pU8DstBufRow;
    }
}

Optimization Points

Reducing duplicate operation
1

uint8_t* pU8SrcBufRowCur = pU8SrcBufRow + i * imgSrc.step;

8. SIMD with reducing naive code

Processing time: 8.2ms

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

__m256i m256iKernelDivTrick = _mm256_set1_epi16(2621);
__m256i m256iZero = _mm256_setzero_si256();
int i32Remain = (imgSrc.cols - 4) % 32;
int i32Back = 32 - i32Remain;

for(int i32Row = 2; i32Row < imgSrc.rows - 2; ++i32Row)
{
    uint8_t* pU8SrcBufRow = pU8SrcBuf + i32Row * imgSrc.step + 2;
    uint8_t* pU8DstBufRow = pU8DstBuf + i32Row * imgDst.step + 2;

    int i32Col = 2;

    for(; i32Col + 32 < imgSrc.cols - 2; i32Col +=32)
    {
        __m256i m256iSumlo = _mm256_setzero_si256();
        __m256i m256iSumhi = _mm256_setzero_si256();

        for(int i = -2; i <= 2; ++i)
        {
            uint8_t* pU8SrcBufRowCur = pU8SrcBufRow + i * imgSrc.step;

            for(int j = -2; j <= 2; ++j)
            {
                __m256i m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRowCur + j));
                __m256i m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero); // latency 1
                __m256i m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero); // latency 1

                m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
                m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);
            }
        }

        m256iSumlo = _mm256_mulhi_epu16(m256iSumlo, m256iKernelDivTrick); // latency 5
        m256iSumhi = _mm256_mulhi_epu16(m256iSumhi, m256iKernelDivTrick);

        __m256i m256iDst = _mm256_packus_epi16(m256iSumlo, m256iSumhi); // latency 3

        _mm256_storeu_si256((__m256i*)pU8DstBufRow, m256iDst);
        pU8SrcBufRow += 32; 
        pU8DstBufRow += 32; 
    }

    pU8SrcBufRow -= i32Back;
    pU8DstBufRow -= i32Back;

    if(i32Back != 0)
    {
        __m256i m256iSumlo = _mm256_setzero_si256();
        __m256i m256iSumhi = _mm256_setzero_si256();

        for(int i = -2; i <= 2; ++i)
        {
            uint8_t* pU8SrcBufRowCur = pU8SrcBufRow + i * imgSrc.step;

            for(int j = -2; j <= 2; ++j)
            {
                __m256i m256iSrc = _mm256_loadu_si256((__m256i*)(pU8SrcBufRowCur + j));
                __m256i m256iSrclo = _mm256_unpacklo_epi8(m256iSrc, m256iZero); // latency 1
                __m256i m256iSrchi = _mm256_unpackhi_epi8(m256iSrc, m256iZero); // latency 1

                m256iSumlo = _mm256_add_epi16(m256iSumlo, m256iSrclo);
                m256iSumhi = _mm256_add_epi16(m256iSumhi, m256iSrchi);
            }
        }

        m256iSumlo = _mm256_mulhi_epu16(m256iSumlo, m256iKernelDivTrick); // latency 5
        m256iSumhi = _mm256_mulhi_epu16(m256iSumhi, m256iKernelDivTrick);

        __m256i m256iDst = _mm256_packus_epi16(m256iSumlo, m256iSumhi); // latency 3

        _mm256_storeu_si256((__m256i*)pU8DstBufRow, m256iDst);
    }
}

Optimization Points

Replace naive code to SIMD

Naive code usually maximizes worst-case, So I should change naive code to simd for average cost if i can

1
2
3
4
5

int i32Remain = (imgSrc.cols - 4) % 32;
int i32Back = 32 - i32Remain;

pU8SrcBufRow -= i32Back;
pU8DstBufRow -= i32Back;
CODE, CODE - SIMD
CODE CODE - SIMD
This post is licensed under CC BY 4.0 by the author.