SIMD 명령어로 두 구의 충돌 검사 최적화 해보기 by object

x86의 SIMD 명령어를 조금 들여다 보다가 재미있는 예가 있어 공유하고자 올린다[1]. 사실 지금까지 SIMD 명령어를 제대로 볼 기회가 없었다. 초보자도 아마 이 글만 제대로 읽으면 어느 정도 감을 잡을 수 있을 것이다.

SIMD는 간단히 말해 하나의 명령어로 여러 데이터를 동시에 처리하는 명령어로 벡터 자료구조를 한번에 계산한다(참고). 예를 들어, 3차원 벡터 두 개를 더 하려면 for 루프를 3번 돌아야 하는데, SIMD 명령어가 있으면 이것을 하나의 명령어로 대체한다. x86은 그 옛날 MMX부터 시작해서 SSE라는 이름으로 SIMD 명령어를 지원한다. 그런데 SSE 버전 4.1부터는 벡터 내적에 관한 명령어가 지원되기 시작했다. 이것으로 충돌 검사 로직 중 일부를 빠르게 할 수 있다.

풀고자 하는 문제는 아래 그림처럼 2차원 공간 상에 두 원이 있을 때, 얘들이 서로 충돌하는지 확인하는 코드이다.

중학교 산수 실력을 빌어 설명하면, 두 원의 충돌 여부는 두 원의 중심 사이의 거리 d가 두 원의 반지름의 합 r1+r2보다 같거나 작으면 충돌이다. 거리 d는 평면 상의 두 원의 중심 A, B를 각각 (x1, y1), (x2, y2)라고 했을 때, 아래처럼 벡터 내적 표현을 빌어 쓸 수 있다. (수식은 간단히 웹 LaTeX를 이용)

따라서 d^2 값과 (r1+r2)^2 값을 비교하면 된다. 이걸 코드로 만들어보자. 단, 여기서 벡터를 2차원이 아니라 3차원으로 확장한다. 3차원 공간에 있는 두 구가 서로 충돌하는지 살펴보는 프로그램이 되겠다. 3차원으로의 확장은 위의 공식에서 z항만 더 추가하면 된다. 원래 [1]에 있는 코드가 약간 불완전해서 수정을 좀 해야 했다.

   1: struct VEC3
   2: {
   3:   float x, y, z;
   4:   float dot() const {
   5:     return x*x + y*y + z*z;
   6:   }
   7: };
   8:  
   9: VEC3 operator-(const VEC3 a, const VEC3 b)
  10: {
  11:   VEC3 result;
  12:   result.x = a.x - b.x;
  13:   result.y = a.y - b.y;
  14:   result.z = a.z - b.z;
  15:   return result;
  16: }
  17:  
  18: struct SPHERE
  19: {
  20:   VEC3  center;
  21:   float radius;
  22: };
  23:  
  24: void Collision_C(const SPHERE& source, const SPHERE targets[], int num,
  25:                  int* collision_results)
  26: {
  27:   for (int i = 0; i < num; ++i)
  28:   {
  29:     float center_distance_sq = (targets[i].center - source.center).dot();
  30:     float radii_sum_sq = (targets[i].radius + source.radius);
  31:     radii_sum_sq = radii_sum_sq * radii_sum_sq;
  32:     if (center_distance_sq < radii_sum_sq)
  33:       ++collision_results[i];
  34:   }
  35: }
  36:  

VEC3 자료 구조를 정의하고 내적(dot product)과 뺄셈을 정의한다. 보다시피 x, y, z에 대해 별도의 계산을 수행한다. 그리고 구(SPHERE) 자료 구조를 만드는데 float 4개가 있는 것과 같다. 중심을 나타내는 3 float 값과 바로 반지름 값이 이어진다. 충돌 검사 함수의 형태는 주어진 구(source)가 targets에 있는 구 배열과 서로 충돌하는지 순회하면서 검사하는 것이다. 코드는 아까 적은 수식을 그대로 옮긴 것 뿐이다.

 

이제 SSE 4.1에서 도입된 백터 내적을 구하는 명령어인 DPPS, DPPD를 이용해서 최적화 해보자. DPPS의 뜻은 Dot Product of Packed Single Precision Floating-Point Values라는 뜻이다. SSE 명령어는 128비트 데이터를 한번에 처리한다. 이 128비트에 float (32비트) 데이터 4개가 포개져서 처리된다. SSE는 4차원 float 벡터에 대한 연산을 한 명령어로 한다. double 값이라면 2차원 벡터에 대한 연산일 것이다. SPHERE 자료 구조는 4개의 float을 가지고 있으므로 SSE 명령을 이용해 최적화 할 수 있다. 아래는 SSE 레지스터(XMMn)의 모습.

일단, 수식부터 다시 정리해보자. 두 구의 자료 구조(중심과 반지름)를 벡터로 표현하면 아래와 같다.

처음에 설명한 것처럼 이 두 점 사이의 거리 d와 반지름의 합을 서로 비교하면 된다. 계산의 편의를 위해 서로 제곱 값을 구하고 그 차이를 구하면 될 것이다. 아래와 같이 k 값을 구해 k 값이 음수라면 두 구는 충돌한 것.

이 k 값을 벡터 a, b로부터 구한다. 그러나 각 항의 부호가 살짝 다르다. x, y, z는 서로 차의 제곱인데 r은 합의 제곱. 그리고 반지름 합의 제곱은 또 앞에 뺄셈이다. 그래서 일부 항의 부호를 살짝 바꾸어서 계산해야 한다. 아래처럼 계산하면 벡터 a, b에서 k를 얻어낼 수 있다.

여기서 a’는 벡터 a에서 마지막 원소의 부호를 뒤집은 것이다. s1은 b에서 a’를 뺀 것이다. 벡터 s2는 s1의 마지막 원소의 부호를 뒤집어 만든다. 그러면 이 s1과 s2를 내적하면 k 값을 얻을 수 있다!!

이 계산 과정을 그대로 SSE 코드로 옮기면 된다. SSE는 바로 기계어로 더 이상 쓸 필요 없이 컴파일러가 제공하는 인트린직을 이용하면 마치 함수처럼 쓸 수 있다. 물론 이 함수는 바로 기계어와 1:1 매핑된다.

   1: #include <smmintrin.h>
   2: void Collision_SSE(const SPHERE& source, const SPHERE targets[], int num,
   3:                    int* collision_results)
   4: {
   5:   __declspec(align(16)) static long _mask[4] = {0, 0, 0, 0x80000000};
   6:   __m128 mask = *(__m128*)_mask;
   7:  
   8:   // a  = {x1, y1, z1,  r1};
   9:   __m128 a  = _mm_load_ps((float*)&source);
  10:   // a1 = {x1, y1, z1, -r1};
  11:   __m128 a1 = _mm_xor_ps(a, mask);
  12:  
  13:   for (int i = 0; i < num; ++i)
  14:   {
  15:     // b = {x2, y2, z2, r2};
  16:     __m128 b = _mm_load_ps((float*)(targets + i));
  17:     // s1 = b - a1 = {x2 - x1, y2 - y1, z2 - z1, r2 - (-r1)};
  18:     __m128 s1 =_mm_sub_ps(b, a1);
  19:     // s2 = {x2 - x1, y2 - y1, z2 - z1, - (r2 + r1)};
  20:     __m128 s2 = _mm_xor_ps(s1, mask);
  21:     // kvec = {(x2 - x1)^2 + (y2 - y1)^2 + (z2 - z1)^2 - (r2 + r1)^2, ... };
  22:     __m128 kvec = _mm_dp_ps(s1, s2, 0xff);
  23:     // k = kvec[0];
  24:     int k = _mm_extract_ps(kvec, 0);
  25:     // 만약 k가 음수라면 최상위 비트가 1이다.
  26:     k >>= 31;
  27:     k &= 1;
  28:     collision_results[i] += k;
  29:   }
  30: }

코드를 설명하면 아래와 같다.

  • 먼저, 라인 9는 SPHERE 자료 구조를 float[4] 형태의 SSE 자료 구조, 128비트 SSE 레지스터 형으로 바꾸는 명령이다. __m128은 128비트 SSE 레지스터를 표현하는 키워드. _mm_load_ps를 이용해 SPHERE의 4개 float을 128비트 SSE 레지스터로 옮긴다(접미사 ps는 packed single).
    • 여기서 주의해야 할 것은 인자 값의 주소가 16바이트에 정렬되어있어야 한다. 즉, source의 주소는 반드시 16의 배수이어야 한다. 5번 라인을 보면 __declspec(align(16))으로 변수가 16바이트에 정렬 맞추어 할당되도록 한다. 정렬이 맞지 않으면 예외가 발생한다.
    • malloc/free 대신 _mm_malloc/_mm_free를 쓰면 된다.
  • 라인 11은 a’를 구한다(코드에서는 a1으로 썼음). 우리가 해야 할 일은 r1의 부호만 바꾸면 된다. 아주 무식하게 한다면 source에 있는 r값의 부호를 바꾸고 라인 9처럼 SSE 레지스터로 바꾸면 될 것이다. 그러나 그렇게 하지 않고 벡터 대 벡터 XOR 연산으로 해결한다.
    • _mm_xor_ps는 이름에서 보듯이 packed single, 즉 float[4] 두 개를 받아 각 원소끼리 XOR 연산을 하는 것. float 데이터 형은 최상위비트가 부호를 나타낸다. 따라서 a에서 마지막 원소 값(반지름)의 최상위비트만 1로 만들어 주면 된다. 그러나 여기서 반지름은 늘 양수이므로 최상위 비트를 반전 시키는 것만 해도 충분할 것이다. 그래서 라인 5의 mask와 벡터 a를 서로 XOR하면 r 값만 부호가 뒤바뀌게 된다. 나머지 x, y, z는 0과 XOR이 되므로 값이 전혀 바뀌지 않는다. (데이터 형이 int라면 +1을 해야 부호가 바뀐다.)
  • 라인 16은 벡터 b를 올리는 코드이다. 역시 배열 targets는 16바이트에 정렬되어야 한다.
  • 라인 18은 s1(=b – a1)를 구하는 코드. _mm_sub_ps는 그 의미를 직관적으로 이해할 수 있다.
  • 라인 20은 s2를 구하는데, 역시 마지막 항의 부호를 바꿔야 한다. 라인 11과 동일한 작업.
  • 라인 22, 24는 이제 s1, s2의 DPPS 명령어로 벡터 내적을 구한다.
    • 0xff라는 마스크 값이 있다. 이것의 자세한 의미는 여기를 보면 나오는데, 특정 원소 값을 계산되지 않게 한다. 이것이 중요한 것은 일종의 분기문 역할을 한다. 벡터 연산 중 조건에 따라 일부 원소의 계산을 막는 것이다. 이 예에서는 큰 의미가 없다.
    • 백터 내적은 사실 스칼라 값이다. 그러나 DPPS는 스칼라 값을 반환하지 않고 내적 값으로 가득 찬 벡터를 반환한다. 여기서 원소 하나를 읽으면 k 값을 최종적으로 얻게 된다. 라인 24가 이 연산을 한다.
  • k 값이 음수인지만 확인하면 충돌 여부를 판가름할 수 있다. 단순하게는 if (k < 0)을 쓸 수 있지만 분기문을 최소화 하는 것은 프로그램 최적화의 기본. k는 비록 int 형으로 선언이 되어있지만 실제로는 float 형이다. 부호 값을 읽는 것은 역시 최상위 비트만 읽으면 되므로 간단한 비트 연산으로 분기문을 대체할 수 있다.

SIMD 명령어를 제대로 써본 적이 없는 나로서는 꽤 흥미로웠다. 인텔의 주장[1]에 따르면, 이렇게 SSE 4.1의 새로운 명령어로 최적화를 하면 C 코드보다 1.5배 빠르다고 한다. 정확하게 C 코드의 기준이 무엇인지 얘는 얼마나 최적화되었는지는 자세히 안 나와있다.

 

SSE는 데이터 정렬이나 배치가 아주 중요하다. 위의 예는 일반적인 4차원 벡터를 선언하고 그것의 배열을 입력 받도록 했다. 이것을 “구조체의 배열”, AOS라고 한다. 그러나 일반적으로 SSE 연산은 이의 반대인 “배열의 구조체”, SOA를 필요로 할 때가 많다. 백터 내적 같은 경우, SSE3 까지만 일반적인 AOS 형태로 주어지면 5개의 SSE 명령이 필요 했다고 한다. 그것을 SSE4에서는 DPPS 하나로 줄였다고 한다.

   1: // Structure of Array (SOA)
   2: struct VECTOR_ARRAY {
   3:   float x[N], y[N], z[N];
   4: } soa;
   5:  
   6: // Array of Strucutre (AOS)
   7: struct VECTOR {
   8:   float x, y, z;
   9: } aos[N];

그런데 컴파일러는 요즘 상당히 똑똑하다. 익히 알다시피 프로그래머가 SSE를 사용하지 않아도 컴파일러가 알아서 SSE 코드로 변환도 해준다. 최적화 옵션만 잘 조절하면 된다. 한번 C 코드를 SSE 명령을 쓰도록 고쳐보았다. 아래 코드를 보면 분명 SSE 코드가 잔뜩 생성됨을 볼 수 있다. 그러나 역시 수학적 원리를 이용한 손수 최적화한 버전에 비하면 못할 수밖에 없다. 아래 기계어 중 굵게 해놓은 것은 한 루프 순환에 해당하는 코드.

Collision_C:
000000013FF610D0  movsxd      r10,r8d 
000000013FF610D3  xor         r8d,r8d 
000000013FF610D6  xor         eax,eax 
000000013FF610D8  test        r10,r10 
000000013FF610DB  jle         Collision_C+6Bh (13FF6113Bh) 
000000013FF610DD  movss       xmm0,dword ptr [rcx+0Ch] 
000000013FF610E2  movss       xmm4,dword ptr [rax+rdx] 
000000013FF610E7  movss       xmm1,dword ptr [rax+rdx+4] 
000000013FF610ED  movss       xmm2,dword ptr [rax+rdx+8] 
000000013FF610F3  movss       xmm3,dword ptr [rax+rdx+0Ch] 
000000013FF610F9  subss       xmm4,dword ptr [rcx] 
000000013FF610FD  subss       xmm1,dword ptr [rcx+4] 
000000013FF61102  subss       xmm2,dword ptr [rcx+8] 
000000013FF61107  mulss       xmm4,xmm4 
000000013FF6110B  mulss       xmm1,xmm1 
000000013FF6110F  addss       xmm3,xmm0 
000000013FF61113  mulss       xmm2,xmm2 
000000013FF61117  addss       xmm4,xmm1 
000000013FF6111B  mulss       xmm3,xmm3 
000000013FF6111F  addss       xmm4,xmm2 
000000013FF61123  comiss      xmm3,xmm4 
000000013FF61126  jbe         Collision_C+5Bh (13FF6112Bh) 
000000013FF61128  inc         dword ptr [r9] 
000000013FF6112B  add         r9,4 
000000013FF6112F  add         rax,10h 
000000013FF61133  inc         r8   
000000013FF61136  cmp         r8,r10 
000000013FF61139  jl          Collision_C+12h (13FF610E2h) 
000000013FF6113B  ret              
000000013FF6113C  nop         dword ptr [rax] 
 
Collision_SSE:
000000013FF61140  movaps      xmm1,xmmword ptr [__xi_z+30h (13FF621A0h)] 
000000013FF61147  movsxd      r8,r8d 
000000013FF6114A  movaps      xmm0,xmmword ptr [rcx] 
000000013FF6114D  xorps       xmm0,xmm1 
000000013FF61150  xor         ecx,ecx 
000000013FF61152  xor         eax,eax 
000000013FF61154  test        r8,r8 
000000013FF61157  jle         Collision_SSE+4Ah (13FF6118Ah) 
000000013FF61159  movaps      xmm3,xmmword ptr [rax+rdx] 
000000013FF6115D  subps       xmm3,xmm0 
000000013FF61160  movaps      xmm2,xmm1 
000000013FF61163  xorps       xmm2,xmm3 
000000013FF61166  dpps        xmm3,xmm2,0FFh 
000000013FF6116C  extractps   r10d,xmm3,0 
000000013FF61173  shr         r10d,1Fh 
000000013FF61177  add         dword ptr [r9],r10d 
000000013FF6117A  add         r9,4 
000000013FF6117E  add         rax,10h 
000000013FF61182  inc         rcx  
000000013FF61185  cmp         rcx,r8 
000000013FF61188  jl          Collision_SSE+19h (13FF61159h) 
000000013FF6118A  ret              
000000013FF6118B  nop         dword ptr [rax+rax] 

참고자료

  1. Improvements in the Intel Core 2 Penryn Processor Family Architecture and Microarchitecture, Intel Technology Journal, Oct 2008 (PDF 내려 받기).

p.s. Windows Live Writer의 Code Snippet으로 코드를 올리니 너무 좋다. RSS에서도 그대로 잘 보인다.


핑백

  • 박피디의 게임 아키텍트 블로그 : Windows Live Writer + Code Snippet 2010-01-17 10:29:32 #

    ... egloos 에서 Code Syntax 보여줄 방법을 계속 못 찾고 있다가, minjang 님 덕분에 좋은 것을 알게 되었다. 바로 Windows Live Writer + Code Snippet 1: void main(){ 2: cout &lt;&l ... more

  • art.oriented : 아이폰 애니메이션 캡처하기 2010-12-22 13:46:32 #

    ... 대해 행렬-벡터 곱셈을 해야 한다. 물론, 이건 가장 무식한 구현이고 실제로 가면 분명 sin/cos 값은 미리 테이블을 만들어 계산할 것이다. SIMD를 이용한 최적화는 사실 기상천외한 수준까지 있으므로 상상 이상의 최적화를 했을 수도 있다. 그러나 이 이야기는 앞서 얘기했듯 비트맵을 가져다가 회전하는 경우에만 해당한다. 위 두 예는 분명 처음 이미지에 대한 ... more

덧글

  • mp3 2010/01/16 16:27 # 삭제 답글

    얼마나 속도가 향상되는지 벤치마킹한 결과는 없나요?
  • object 2010/01/16 16:30 #

    >> 이렇게 SSE 4.1의 새로운 명령어로 최적화를 하면 C 코드보다 1.5배 빠르다고 한다.

    참고자료 [1]에 살짝 비교 되어있는데, 일단 코드만 보더라도 분명 SSE4.1이 빠르겠죠. 그런데 '얼마나' 더 빠른가는 이게 워낙 가변적인 요소가 많아서 직접 해보셔야 할 겁니다. 일단 눈으로봐도 1.5배 정도 차이가 나기는 나겠습니다.
  • Ya펭귄 2010/01/16 16:30 # 답글

    Matlab m-script는 아예 스크립트 차원에서 SIMD를 지원해주던 기억이 나네요.... (애초에 Matlab의 자료형이 무조건 double-float이니 당연하다면 당연하겠지만....) 차량의 구동륜에 대한 마찰력 해석이었는데 그게 알고리즘상으로 바퀴 네 개에 대한 병렬연산이 가능했던지라 일반 스칼라 연산 코드를 4-order벡터연산 코드로 바꾸는 순간 실행속도가 2.5배 정도 증가하더라는....(아직 Matlab이 내부적으로 double-float를 x87 FPU에서만 연산하던 시절에도 속도증가가 있더군요...)

    그런데 예전에 어디에서인가에서 x87 FPU의 double-float가 80비트 확장정밀도인지라 연산시 64bit float 가 overflow되는 경우 overflow되는 부분에 대한 손실 없이 처리하도록 만드는 게 가능했다가 SSE2가 도입되면서 그런 식으 처리가 좀 어려워졌다는 이야기를 읽은 듯 합니다.
  • object 2010/01/16 16:36 #

    데이터 패러럴한 작업은 SIMD로 크게 성능이 나아지겠죠. 2.5배면 상당한 수준이군요.

    x87은 말씀대로 80비트, 아이태니엄도 64비트 플로팅 레지스터가 정밀도 때문에 82비트로 확장해서 가지고 있는데요. SSE의 xmm 레지스터나 AVX의 ymm은 그러한 이야기가 잘 안보이네요. 저도 그런 계산 쪽은 잘 아는 바가 없어서;;
  • DJ™ 2011/06/06 09:39 # 삭제

    SSE의 경우에는 x87처럼 정밀도를 위한 추가장치가 없는것으로 알고 있습니다.
    그래서 x87에서의 연산 결과와 SSE의 연산 결과가 달라지는 경우도 종종 있더군요..
  • 최종욱 2010/01/16 19:52 # 답글

    푸하하... 앞으로는 더 미세한 부분의 알고리즘들도 SIMD, GPU 등에 맞추어 다시 쓰여질 것 같습니다. 재미있네요. 잘 봤습니다!
  • object 2010/01/17 04:31 #

    우리가 몰라서 그렇지 이미 많이들 쓰고 있죠. SSE 4.1에 이런 명령어가 들어온 것도 필요가 있어서였습니다. 비슷하게 SSE 4.2부터는 텍스트 파싱을 편하게 해주는 함수도 추가되었죠.
  • ... 2010/01/17 04:21 # 삭제 답글

    잘 봤습니다.
    성능향상이 약 1.5배 라고 쓴걸 보니 궁금한 점이 생기네요

    SSE로 최적화 했는데 성능향상이 굉장히 크다면
    그 알고리즘은 좋은 알고리즘이라고 판단 해야 할까요?
    아니면 나쁜 알고리즘이라고 판단해야 할까요??
  • object 2010/01/17 04:30 #

    알고리즘 중에는 병렬성을 찾기 쉬운 것도 있고 이 예도 그러한 것일 뿐입니다. 좋고 나쁘고를 이야기 할 것은 아닙니다. SSE로 이득을 볼 수 있는 것은 벡터 타입 같은 데이터 패러럴한 작업이 가능할 때입니다. GPGPU도 마찬가지이고요.
  • summerligh 2010/01/18 02:14 # 삭제 답글

    FFT 같은 경우는 SSE 명령어를 잘 쓰면 3배 가까이도 빨라지더군요. 재귀의 80%가 죄다 곱셈-덧셈의 반복이라...
  • 박PD 2010/01/18 10:02 # 답글

    덕분에 좋은 거 배웠습니다. 감사합니다.
  • Ciel 2010/01/19 15:36 # 삭제 답글

    너무 유익한 내용이네요.
  • 미친감자 2010/01/25 17:29 # 답글

    지송합니다. 항상...읽다가 어려워서 ㅡ,.ㅡ;;;; 말아버리내요...
    다시 들어와서 차분하게 천천히 잘 읽어보겠습니다...
  • 실버테란 2010/12/21 08:40 # 답글

    안녕하세요. 블로그를 보고 Code Snippet을 사용하려고 하는데요. 이글루스에 사용하면 자동으로 <br>이 들어가는 현상이 발생하지 않으신가요? 전 지금 그거때문에 수동으로 제거해주고 있어서요. ㅠㅠ
    좋은 게시물 잘 보고 갑니다.
  • yanus 2011/01/19 18:24 # 삭제 답글

    FFT관련 정보를 찾다가 우연히 SIMD에 관한글을 읽고 들어왔는데, 좋은 글 잘 읽고갑니다.

    SIMD고속화에 대해 한가지 팁이라고 하면,
    때로는 SIMD를 적용하기 위해서 SIMD에 맞게 바꾸는 경우도 있습니다. 물론 같은 동작을 하게 하는것이지만, SIMD로 가장 잘 동작하게 만들게 하는것이죠. (메모리접근성, 파이프라인 깨지지 않도록 하기 위해)
    실제로 파이프라인을 깨지지 않게 하기 위해서, 명령어를 늘리는 경우도 있습니다.

    그리고 일반적으로 수치연산작업만을 SIMD 처리하는 것에 많은 촛점을 두는 경향이 있지만 실은 알고리즘에 의해 꽤 여러가지 것들을 SIMD를 사용하여 처리할 수 있습니다. (그런 알고리즘들이 보면 이상한 알고리즘이 꽤 많지요... 하지만, 고속화에 초점을 두다 보면 그런것들을 만들어야지 어쩔수 없답니다... OTL)

    좋은 글 잘 보고 갑니다.~ 좋은 하루 되세요~
댓글 입력 영역