여기서 SIMD 명령어를 사용하면 한 번에 더 많은 요소를 처리하는 데 어떻게 도움이 될지 생각해 보십시오. x86 환경이라면 Intel Intrinsics Guide에서 SSE 계열 내장 함수를 살펴보고, 어떤 것이 이 경우에 적합할지 판단해 보십시오. 현재 코드는 8비트 정수를 한 번에 하나씩 다룬다는 점을 떠올려 보십시오.

이 코드는 슬라이딩 윈도를 사용해 범위 [input.data() + pos - radius, input.data() + pos + radius]에 걸쳐 위치 pos의 블러를 계산하며, 너비가 2 * radius를 넘지 않도록 슬라이딩 윈도에서 input.data() + pos + radius - 1를 제거합니다. _mm_loadu_si64 SSE 명령어로 64비트(즉, 8개의 pos 시작점)를 한 번에 로드해 보십시오.

스칼라 방식의 문제는 누적 합(prefix sum)을 사용합니다. 이를 SIMD 방식으로 처리하려면 벡터화된 누적 합 알고리즘도 알아야 합니다. 벡터 [1, 2, 3, 4]의 경우 다음과 같습니다.

1. 기준값에서 시작합니다: res = [1, 2, 3, 4]
2. 기준값을 오른쪽으로 2⁰ = 1만큼 시프트하고 결과에 더합니다: res = [1, 2, 3, 4] + [0, 1, 2, 3] = [1, 3, 5, 7]
3. 기준값을 오른쪽으로 2¹ = 2만큼 시프트하고 더합니다: res = [1, 3, 5, 7] + [0, 0, 1, 2] = [1, 3, 6, 9]
4. 기준값을 오른쪽으로 2² = 4만큼 시프트하고 더합니다: res = [1, 3, 6, 9] + [0, 0, 0, 0] = [1, 3, 6, 9]
5. 더 이상 수행할 덧셈이 없으므로 멈춥니다(덧셈이 더 있다면 2³ = 8만큼 시프트하는 식으로 계속합니다).

이 내장 함수들은 리틀 엔디언 이진수를 다루므로, 내장 함수로 수행하는 시프트 연산은 오른쪽이 아니라 왼쪽 시프트여야 한다는 점에 유의하십시오(위의 알고리즘은 빅 엔디언 표현을 사용합니다). 이 왼쪽 시프트에는 _mm_slli_si128를 사용하게 됩니다.

8개 단위 블록에 들어맞지 않아 SIMD 명령어로 읽을 수 없는 나머지 요소는 어떻게 처리해야 하겠습니까?

배경:

이 실습의 이미지 스무딩 알고리즘은 슬라이딩 윈도 방식을 사용해, 어떤 기준점을 중심으로 한 반경 안에 있는 값들의 누적 합을 계산합니다. 목표는 스칼라 누적 합을 SIMD 누적 합으로 바꿔 계산을 훨씬 빠르게 수행하는 것입니다.

벡터 양의 누적 합을 계산하려면 먼저 벡터에 있는 값들에서 시작합니다. 그런 다음 벡터를 오른쪽으로 한 요소만큼 이동시킨 후 누적합니다. 두 번째 반복에서는 2만큼 시프트하고, 그다음 반복에서는 4만큼 시프트하는 식으로 진행합니다. 따라서 이 알고리즘은 순차적 누적 합에서 볼 수 있는 선형 복잡도가 아니라 로그 복잡도를 가집니다.

컴파일러가 이를 위한 SIMD 코드를 생성하기는 상당히 어려우므로, 같은 효과를 얻으려면 우리가 직접 컴파일러 내장 함수를 사용해야 합니다.

먼저 원래 코드의 벤치마크를 측정해 봅시다:

------------------------------------------------------------------------------
Benchmark                                    Time             CPU   Iterations
------------------------------------------------------------------------------
bench_partial_sum/iterations:100000       25.4 us         25.3 us       100000

해결책

SSE(Streaming SIMD Extensions) 풀이부터 시작해 봅시다. 슬라이딩 윈도 누적 합에서 출력 배열의 pos 위치에 도달하면, const uint8_t subtractVal = input.data() + pos - radius - 1이 가리키는 위치의 값(윈도 범위에서 빠져나간 값)을 합계에서 빼야 합니다. 그런 다음 슬라이딩 윈도우에 방금 추가된 새 값인 const uint8_t addVal = input.data() + pos + radius를 더해야 합니다.

먼저, 크기가 [1, 2 * radius + 1]인 슬라이딩 윈도에서 누적 합을 쌓아 얻은 currentSum을 SIMD 연산을 수행할 수 있는 형태로 변환해야 합니다.

_mm_set1_epi16을 사용하여 16비트 currentSum의 복사본 여덟 개를 담은 128비트 벡터 값을 만듭니다. 이를 입력을 브로드캐스팅한다고 합니다.

__m128i current = _mm_set1_epi16(currentSum); // create vector of eight copies of currentSum

그런 다음 반복마다 데이터의 요소 여덟 개를 처리합니다. 필요한 값들은 다음과 같이 변환합니다:

• subtractVal - _mm_loadu_si64를 사용하여 입력에서 다음 uint8_t 여덟 개를 128비트 벡터의 처음 64비트에 로드해야 합니다.
• addVal - subtractVal과 같은 내장 함수(_mm_loadu_si64)를 사용해야 합니다.

이 값들을 벡터 값에 로드하고 나면, 1바이트 정수 여덟 개가 이 벡터들의 처음 64비트에만 채워져 있다는 문제가 생깁니다. 위에서 만든 current 벡터와 이 벡터들 간에 산술 연산을 수행할 수 있도록, _mm_cvtepu8_epi16 내장 함수를 사용하여 이 1바이트 정수들을 2바이트 정수로 확장해야 합니다.

지금까지의 코드는 다음과 같습니다:

const uint8_t* subtract_ptr = input.data() + pos - radius - 1;
const uint8_t* add_ptr = input.data() + pos + radius;
const uint16_t* output_ptr = output.data() + pos;

__m128i current = _mm_set1_epi16(currentSum);

int i = 0;
for (; i + 7 + pos < limit; i += 8) {
    __m128i sub_packed_8 = _mm_loadu_si64(subtract_ptr + i); // load the eight `uint8_t` integers in range [subtract_ptr, subtract_ptr + 8)
    __m128i add_packed_8 = _mm_loadu_si64(add_ptr + i); // load the eight `uint8_t` integers in range [add_ptr, add_ptr + 8)
    __m128i sub_ext_16 = _mm_cvtepu8_epi16(sub_packed_8); 
    __m128i add_ext_16 = _mm_cvtepu8_epi16(add_packed_8);
}

슬라이딩 윈도우를 앞으로 이동시킬 때, 여덟 요소 각각에 해당하는 addVal을 더하고 subtractVal을 뺍니다. 이것이 차이값(diff)이 됩니다. 따라서 이 차이에 해당하는 SIMD 연산은 _mm_sub_epi16(add_ext_16, sub_ext_16)입니다.

지금부터는 __m128i 타입을 나타내기 위해 v8i를 사용하겠습니다(v8i는 "요소 8개짜리 벡터"를 의미합니다):

using v8i = __m128i;
...
v8i diff = _mm_sub_epi16(add_ext_16, sub_ext_16);

그런 다음 벡터 누적 합 알고리즘을 수행합니다. 한 가지 중요한 점에 주의하십시오! 이 내장 함수들은 리틀 엔디언 이진수를 제공하는데, 이는 최상위 바이트가 가장 큰 메모리 주소에 위치한다는 뜻입니다. 앞서 설명한 알고리즘에서는 오른쪽 비트 시프트를 한다고 했지만, 이는 빅 엔디언 표현에 해당하는 이야기였습니다. 예를 들어, 32는 빅 엔디언 8비트 표현으로 10000000이지만, 리틀 엔디언으로는 00000001입니다.

따라서 이 알고리즘에서는 처리 중인 16비트 정수 8개에 걸쳐 모든 diff 기여분의 누적 합을 쌓을 때 왼쪽 시프트를 수행해야 합니다.

// Copied from above
v8i diff = _mm_sub_epi16(add_ext_16, sub_ext_16);

// Calculate the prefix sum of all the `diff`s for the eight elements:
v8i delta = _mm_add_epi16(diff, _mm_slli_si128(diff, 2)); // shift `diff` by 2 bytes to the left (i.e. drop one 16-bit integer) and add to `diff`
delta = _mm_add_epi16(delta, _mm_slli_si128(delta, 4)); // shift `delta` by 4 bytes to the left (i.e. drop two 16-bit integers) and add to `delta`
delta = _mm_add_epi16(delta, _mm_slli_si128(delta, 8)); // shift `delta` by 8 bytes to the left (i.e. drop four 16-bit integers) and add to `delta`

그런 다음 이 delta를 기준 버전 current에 합쳐 i번째 반복의 result를 얻습니다:

v8i result = _mm_add_epi16(delta, current);
_mm_storeu_si128((v8i*)(output_ptr + i), result); // store result into memory at address *(output_ptr + i)

마지막으로, 16비트 정수 여덟 개 중 마지막 값을 가져와 다음 루프 반복을 위해 브로드캐스트합니다:

current_sum = (uint16_t)_mm_extract_epi16(result, 7); // extract last element
current = _mm_set1_epi16(current_sum); // broadcast current_sum into 128-bit register

루프가 끝나면 (출력을 쓰기 위한) pos를 모든 반복에서 누적된 i 값만큼 증가시켜야 합니다:

pos += i;

지금까지의 내용을 모두 종합하면 다음과 같습니다:

  // --- THE CODE ABOVE HERE IS UNCHANGED ---
  ...
  // 2. main loop.
  limit = size - radius;
  const uint8_t* subtract_ptr = input.data() + pos - radius - 1;
  const uint8_t* add_ptr = input.data() + pos + radius;
  const uint16_t* output_ptr = output.data() + pos;

  using v8i = __m128i;
  v8i current = _mm_set1_epi16(currentSum);

  int i = 0;
  for (; i + 7 + pos < limit; i += 8) {
    v8i sub_pckd_8 = _mm_loadu_si64(subtract_ptr + i);
    v8i add_pckd_8 = _mm_loadu_si64(add_ptr + i);
    v8i sub_ext_16 = _mm_cvtepu8_epi16(sub_pckd_8);
    v8i add_ext_16 = _mm_cvtepu8_epi16(add_pckd_8);

    v8i diff = _mm_sub_epi16(add_ext_16, sub_ext_16);

    v8i delta = _mm_add_epi16(diff, _mm_slli_si128(diff, 2));
    delta = _mm_add_epi16(delta, _mm_slli_si128(delta, 4));
    delta = _mm_add_epi16(delta, _mm_slli_si128(delta, 8));

    v8i result = _mm_add_epi16(delta, current);
    _mm_storeu_si128((v8i*)(output_ptr + i), result);

    currentSum = static_cast<uint16_t>(_mm_extract_epi16(result, 7));
    current = _mm_set1_epi16(currentSum);
  }
  pos += i;
  ...
  // --- THE CODE BELOW HERE IS UNCHANGED ---

이러한 x86 컴파일러 내장 함수를 사용하면 벤치마크에서 코드가 두 배 이상 빨라집니다!

------------------------------------------------------------------------------
Benchmark                                    Time             CPU   Iterations
------------------------------------------------------------------------------
bench_partial_sum/iterations:100000       10.4 us         10.4 us       100000

해답 링크