never knows best

한번 SIMD를 사용해서 행렬곱셈을 해 봤습니다.
있는 코드가 있으면 편하게 긁어다 붙이려 했는데 구글을 아무리 뒤져도 SIMD를 이용해서 행렬*벡터만 보이고 행렬*행렬을 하는 코드가 안 보이길래 만들어버렸군요. 혹시 영어로 검색하는 사람들을 위해서 영어로도 써 봅니다.
Optimized matrix multiplication using SIMD
정도면 되려나요.


* SIMD

여러개의 계산을 병렬처리하는 방법으로 아마 펜티엄 2인가 3에서 128비트 레지스터에 부동소수점 4개를 올려서 한방에 계산하는 뭐 그런게 나왔던것 같습니다. 서로 의존성이 없는 값들을 계산하기에 적합한데 특히 그래픽스에서 행렬과 벡터를 다룰 때 유용하게 사용할수 있습니다.


* 아이디어

기본적으로 SIMD를 이용하면 실수 4개의 덧셈 혹은 곱셈을 한번에 할 수 있습니다.
행렬 계산식을 잘 보면 다음의 두 행렬을 곱한다고 쳤을 때
( a11 a12 a13 a14 )   ( b11 b12 b13 b14 )
( a21 a22 a23 a24 )   ( b21 b22 b23 b24 )
( a31 a32 a33 a34 )   ( b31 b32 b33 b34 )
( a41 a42 a43 a44 )   ( b41 b42 b43 b44 )

두 행렬 A, B를 곱해서 나오는 행렬 C의 첫행은

c11 = a11 * b11 + a12 * b21 + a13 * b31 + a14 * b41
c12 = a11 * b12 + a12 * b22 + a13 * b32 + a14 * b42
c13 = a11 * b13 + a12 * b23 + a13 * b33 + a14 * b43
c14 = a11 * b14 + a12 * b24 + a13 * b34 + a14 * b44

이므로

중간중간에 더하기를 기준으로 잘라보면

(c11, c12, c13, c14) =
(a11 * (b11, b12, b13, b14)) +  
(a12 * (b21, b22, b23, b24)) +  
(a13 * (b31, b32, b33, b34)) +  
(a14 * (b41, b42, b43, b44))

가 됩니다.

여기서 SIMD를 사용하면 a11 * (b11, b12, b13, b14) 를 한 명령어로 처리할 수 있으므로 속도 증가를 기대할 수 있습니다.



* VC 라이브러리

VC에서는 SIMD를 위해서 쓰기 편하게 랩핑된 클래스를 제공합니다. 인라인 어셈이나 알아보기 불편한 함수 셋트가 있긴 한데 편한걸 두고 굳이 힘들게 쓸 이유는 없겠죠;
VC에서 128비트 레지스터에 올라갈 변수는 __m128형입니다. 단순히 float[4]라고 보시면 되겠네요. 한편 SIMD계산을 위한 클래스는 F32vec4입니다. F32vec4는 <fvec.h>를 인클루드해야만 사용가능합니다.



* 코드

struct CMatrix
{
    union
    {
       struct 
       {
            __m128 _L1, _L2, _L3, _L4;
       };
       struct 
       {
            float _11, _12, _13, _14;
            float _21, _22, _23, _24;
            float _31, _32, _33, _34;
            float _41, _42, _43, _44;
        };
    };
    void    MultiplyWith(const CMatrix& mat)
    {
        F32vec4 t1, t2, t3, t4;

        t1  = mat._L1 * F32vec4(_11);
        t1 += mat._L2 * F32vec4(_12);
        t1 += mat._L3 * F32vec4(_13);
        t1 += mat._L4 * F32vec4(_14);

        t2  = mat._L1 * F32vec4(_21);
        t2 += mat._L2 * F32vec4(_22);
        t2 += mat._L3 * F32vec4(_23);
        t2 += mat._L4 * F32vec4(_24);

        t3  = mat._L1 * F32vec4(_31);
        t3 += mat._L2 * F32vec4(_32);
        t3 += mat._L3 * F32vec4(_33);
        t3 += mat._L4 * F32vec4(_34);

        t4  = mat._L1 * F32vec4(_41);
        t4 += mat._L2 * F32vec4(_42);
        t4 += mat._L3 * F32vec4(_43);
        t4 += mat._L4 * F32vec4(_44);

        _L1 = t1; _L2 = t2; _L3 = t3; _L4 = t4;
    }
};




* 주의사항

__m128의 경우 반드시 16byte로 정렬된 메모리에 올라가 있어야만 합니다. 내부적으로 __declspec(align(16)) 을 사용하는데, 이렇게 정렬된 메모리상에 있는 변수의 경우 함수의 인자로 넘길때 값에 의한 (call-by-value) 전달은 허용되지 않습니다. 보통은 const 레퍼런스를 넘기를 방법으로 회피할 수 있죠.
또한 std::vector를 사용할 때 애로사항이 꽃피는데 VC의 std::vector 구현 중 resize() 내부에서 신나게 값에 의한 전달을 하기 때문입니다. 이 현상을 가장 간단하게 회피하는 방법은 std::vector를 대신하는 전용 컨테이너를 만드는것이죠. 이때 주의할 점은 단순히 new로 CMatrix를 할당했을 경우 정렬된 메모리를 얻을 수 없기 때문에 _L1에 엑세스하려는 순간 access violation error가 뜨면서 코드가 죽는 문제가 생깁니다.
이 문제를 회피하기 위해서는 new, new[], delete, delete[]를 모두 오버로드해서 내부적으로 _aligned_malloc()과 _aligned_free()를 사용하도록 신경써야 합니다. 또한 CMatrix를 포함하고 있는 클래스가 있을 경우 그 또한 정렬된 메모리상에 올라가야 합니다. 당연한 얘기죠?



* 퍼포먼스

제 테스트로는 제 CPU(P4 E6300 듀오) 에서 단순히 곱셈/덧셈으로 구현한 행렬곱보다 80-90%정도의 성능향상이 있었습니다.



* 참고자료

http://www.gamasutra.com/features/20000131/barad_03.htm
http://www.gamedev.net/community/forums/topic.asp?topic_id=526505