[Parallel Computing] (16) Tiling for Matrix Multiply

Strip Mining은 한 개의 Loop를 여러개로 나눠서 수행하는 것을 말한다. 보통 Vector 연산을 수행하기 위해서 사용한다.

Matrix Mulitply에서도 Strip Mining을 적용해서 사용할 수 있다. 아래와 같은 Matrix 연산을 생각해보자.  

for (int i = 0; i < N; i++)
  for (int j = 0; j < N; j++)
    for (int k = 0; k < N; k++)
        c[i][j] += a[i][k] + b[k][j];

'a' 행렬과 'b' 행렬를 곱하여 'c' 행렬에 저장하는 동작이다.
여기서 잘 살펴봐야하는 건, "'a'와 'b' 행렬이 어떻게 메모리에 저장되어 있는 가" 하는 점이다. 코드상에는 2D 배열로 표현되지만 실제 저장은 1D형태로 저장된다. 일반적으로 2D 배열은 하나의 전체 Row를 순차적으로 저장하고, 다음 Row를 저장하는 형태로 메모리에 저장된다. 즉, 한번의 메모리 접근으로 동일한 Row의 여러 Column들을 읽어 올 수 있다. 'a' 행렬의 경우, k for-loop이 수행될 때마다 같은 Row의 여러 Column을 접근하므로 한번의 Cache 접근(1 Cache Miss)으로 여러 번의 Cache Hit를 발생시킬 수 있다. 반대로 'b' 행렬의 경우, k for-loop이 수행될 때마다 다른 Row의 같은 Column을 접근하므로 한번의 Cache 접근(1 Cache Miss)으로 한번의 데이터만 읽어올 수 있다.
물론, 'b' 행렬에 접근할 때, 이전 k for-loop에서 접근한 데이터들이 Cache에 남아 있어, 다음 k for-loop에서 같은 Row의 다른 Column을 접근하면 Cache Hit이 발생할 수 도 있다. 이는 'b' 행렬 전체를 Cache에 저장할 수 있어야 한다는 말과 같다. 'b' 행렬이 충분히 크게되면 Cache 크기를 넘어갈 수 있다.
이를 수치로 표현해보기 위해서 Cache 내에는 'a' 행렬를 저장하기 위한 1개의 Cache line과 'b' 행렬를 저장하기 위한 1개 Cacheline이 있다고 가정해보자. 이 때 Cacheline의 크기(m)은 행렬의 크기(N) 보다 작다고 가정한다.
'a'행렬의 Cache miss는 Cache에 없을 때 한번 발생하고, 나머지 m-1에 대해서는 발생하지 않는다. 즉, m번마다 한번 발생한다. 총 'a'행렬의 접근 횟수는 N*N이므로
'a' 행렬의 Cache miss는 N(i for-loop)*N(k for-loop)/m 으로 생각할 수 있다.
'b'행렬의 Cache miss는 Cache에 없을 때마다 발생한다. 서로 다른 Row의 같은 Column은 for-loop를 수행할 때마다 Cache에 존재하지 않으므로 'b' 행렬를 접근할 때마다 발생한다.
'b' 행렬의 Cache miss는 N(i for-loop)*N(j for-loop)*N(k for-loop)이다.
따라서 총 Cache miss는 N^2/m + N^3 이다.

이를 더 효율적으로 할 수 있는 방법은 없을 까?
그에 대한 해답은 위에서 언급한 것에서 힌트를 찾아볼 수 있다.

물론, 'b' 행렬에 접근할 때, 이전 k for-loop에서 접근한 데이터들이 Cache에 남아 있어, 다음 k for-loop에서 같은 Row의 다른 Column을 접근하면 Cache Hit이 발생할 수 도 있다. 이는 'b' 행렬 전체를 Cache에 저장할 수 있어야 한다는 말과 같다. 'b' 행렬이 충분히 크게되면 Cache 크기를 넘어갈 수 있다.

즉, 'b' 행렬이 Cache에 들어올 수 있도록 for-loop의 연산을 나누는 것이다. 이를 tiling 기법이라고 한다. tiling 기법을 적용하면 아래와 같은 코드로 나타낼 수 있다.

// iterations of subblocks
for (int ii = 0; ii < N; ii += B)
  for (int jj = 0; jj < N; jj += B)
    for (int kk = 0; kk < N; kk += B)
       // B x B submatrix multiplications
       for (int i = ii; i < MIN(ii+B, N); i++)
          for (int j = jj; j < MIN(jj+B, N); j++)
             for (int k = kk; k < MIN(kk+B, N); k++)
                 c[i][j] += a[i][k] * b[k][j];

아래 마지막 3개의 for-loop은 subblock multiplication으로 작은 단위의 Matrix Multiply를 수행한다. 이 때, 하나의 subblock이 Cache에 모두 저장될 수 있다고 생각해보자.
'a', 'b' subblock 모두 Cache에 있으므로 접근 횟수를 Cacheline 크기로 나눈 값이 Cachemiss 발생횟수이다.
즉, 각각 Cachemiss는 B^2/m이므로 총 Cachemiss는 2*B^2/m으로 표현될 수 있다.
이러한 Cachemiss는 가장 바깥의 3개의 for-loop 수행만큼 발생한다.
1개의 for-loop이 N/B 만큼 수행하므로 (N/B)^3만큼 수행된다.
두 값을 곱하게 되면 총 Cachemiss를 알 수 있다.
2*B^2/m * (N/B)^3 = 2N^3/(B*m) 으로 계산된다.

subblock으로 나누지 않았을 때와 비교하면 그 차이가 보인다.
N^2/m + N^3 vs. 2*N^3/(B*m)

하나의 Cache에 담을 수 있는 Blocksize가 클 수록 Cachemiss 횟수가 줄어듬을 알 수 있다.