공간 지역성을 높이기 위한 루프 재배치

행렬의 곱셈을 구하는 함수는 대개 세 개의 중첩된 루프를 사용하여 구현하며, 이들은 각각 인덱스 i, j, k로 식별된다. 이 루프들의 순서를 바꾸고 일부 코드를 수정하면, 하나의 동일한 행렬 곱을 수행하는 코드라고 하더라도 여러가지 코드 형태를 보일 수 있다. 이때 각각의 버전은 루프 순서에 의해 구분할 수 있다. 해당 글에서는 총 6가지의 버전을 비교하며 작성되었다.

상위 레벨에서는 어떤 버전으로 코드가 작성되었든 비슷하다. 만일 덧셈의 결합성이 성립한다면 각 버전들은 동일한 결과를 낼 것이다. (부동소수점 덧셈은 교환성은 성립되지만, 결합성이 성립되지 않는다는 점을 기억하자) 각 버전은 O(n^3)의 전체 연산을 수행하며, 동일한 수의 덧셈과 곱셈을 수행한다. A와 B의 n^2 개의 원소들은 n번 읽히게 될 것이며, C의 n^2개의 원소 각각은 n개의 값들을 더해서 계산될 것이다. 그러나 가장 안쪽 루프의 동작을 분석할 경우 접근 횟수와 지역성에 있어서 차이가 있음을 발견할 수 있다. 이 해석을 위해서 다음의 가정을 한다.

• 각 배열은 double의 n x n 배열이며, sizeof(double) == 8 이다.
• 32바이트 블록 크기를 갖는 (B = 32) 단일 캐시를 가지고 있다.
• 배열 크기 n은 매우 커서 단일 행렬의 행은 L1 캐시에 들어가지 않는다.
• 컴파일러는 지역변수를 레지스터에 저장하고, 따라서 루프 내 지역변수에 대한 참조는 load나 store 연산을 사용하지 않는다(지역 변수의 load에는 메모리의 개입이 없다).

들어가기에 앞서서

C의 array는 row-major하게 할당되어 있다.

때문에 column base로 접근하느냐 row base로 접근하느냐에 따라 Miss rate가 극단적으로 변화한다.

먼저 column이 변화하고 한 줄의 column을 다 읽었을때 새로운 row로 루프가 진행되는 경우를 살펴보자

for (i = 0; i < n; i++)
	sum += a[0][i];

이 경우에는 연속적인 원소에 접근하게 된다. 이 경우에 Miss rate가 가장 커지려면 블록의 사이즈가 원소의 사이즈보다 커서 공간 지역성을 활용할 수 있다.  이 경우의 miss rate는 다음과 같다. 

전체 B 바이트를 읽을 때마다 (단일 매핑 캐시의 전체), sizeof(aij) 만큼의 캐시 미스가 발생하였다. 이때 sizeof(aij)는 모든 원소가 같은 크기라고 가정하였다.

반면, row를 가장 깊은 루프에서 변화시키며 한줄의 row를 다 읽었을 때 다음과 같이 새로운 column으로 루프가 진행될 경우, C의 배열 설계와 정반대로 가서 동떨어져 있는 배열 원소에 계속해서 접근하게 된다. 

for (i = 0; i < n; i++)
	sum += a[i][0]

그리하여 공간 지역성이 성립이 되지 않아 Miss rate가 1이 되어 캐시를 사용하지 않는 것과 같은 역할을 수행한다.

# Case 1 (i- j-k)

for (i = 0; i < n; i++)
{
    for (j = 0; j < n; j++) {
    	sum = 0.0
        for (k = 0; k < n; k++)
        	sum += a[i][k] * b[k][j];
        c[i][j] = sum;
    }
}

해당 케이스는 클래스 AB로 분류된다. (가장 안쪽의 루프에서 [i][k] -> [k][j] 로 이어지기 때문)

이때, a는 stride 1으로서 sizeof(a) / B의 miss rate 를 갖게 되고, 앞서 가정하였듯 블록의 크기가 32 바이트이므로 각각의 원소가 double(8 bytes)라는 점에서 0.25의 miss rate를 갖게된다.

b의 경우에는 stride n으로 스캔이 이뤄지기 때문에, 1의 miss rate를 갖는다.

총 1.25 번의 미스가 도출되었다.

# Case 2 (j- i-k)

for (j = 0; j < n; i++)
{
    for (i = 0; i < n; j++) {
    	sum = 0.0
        for (k = 0; k < n; k++)
        	sum += a[i][k] * b[k][j];
        c[i][j] = sum;
    }
}

해당 코드는 직전의 Case #1과 동일한 클래스 AB로 분류된다. 앞서 a가 0.25의 miss rate를 갖던 것과 마찬가지로  stride 1으로 스캔이 이뤄져 0.25의 miss rate를 갖는다.

b의 경우도 앞의 이유에서처럼 stride n을 가지므로 1의 miss rate를 갖는다.

총 1.25 번의 미스가 도출되었다.

# Case 3 (j- k-i) 와 # Case 4 (k-j-i)

for (j = 0; j < n; j++)
{
    for (k = 0; k < n; k++) {
        r = B[k][j];
        for (i = 0; i < m; i++)
            C[i][j] += A[i][k] * r;
    }
}

for (k = 0; k < n; j++)
{
    for (j = 0; j < n; k++) {
        r = B[k][j];
        for (i = 0; i < m; i++)
            C[i][j] += A[i][k] * r;
    }
}

클래스 AC이다.

다음의 연산에는 두가지 문제가 있다. 앞서 클래스 AB가 load 연산만을 2번 사용했던 것과 달리 해당 코드에서는 store 연산 2번, load 연산 2번을 사용하고 있다. 

또한 내부 루프에서 A와 C의 열들을  stride n으로 실행한다(가장 깊은 루프에서 i가 바뀌면서 루프 진행). 그 결과는 매 load마다(A와 C) 한 번의 미스(miss rate = 1)가 발생하여 총 2번의 미스가 매 실행마다 발생한다.                                                              

# Case 5 (i- k-j) 와 Case 6 (k-i-j)

for (i = 0; i < m; i++) {
    for (k = 0; k < n; k++) {
        r = A[i][k];
        for (j = 0; j < n; j++)
            C[i][j] += r * B[k][j];
    }
}

for (k = 0; k < n; k++)
{
    for (i = 0; i < n; i++) {
        r = A[k][i];
        for (j = 0; j < n; j++)
            C[i][j] = r * B[k][j];
    }
}

클래스 BC이다.

직전의 클래스 AC와 동일하게 2번의 load, store 연산이 각각 수행되고 있으나, 가장 깊은 루프에서 C와 B가 행우선으로 stride-1 접근 패턴을 적용하기 때문에 각각의 load 연산에서 0.25의 miss rate를 갖는다.

총 0.5의 miss rate를 갖는다.