행렬 곱셈 : 행렬 크기의 작은 차이, 타이밍의 큰 차이

행렬 곱셈 : 행렬 크기의 작은 차이, 타이밍의 큰 차이

---

다음과 같은 행렬 곱셈 코드가 있습니다.

for(i = 0; i < dimension; i++)
    for(j = 0; j < dimension; j++)
        for(k = 0; k < dimension; k++)
            C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];

여기서 행렬의 크기는로 표시됩니다 dimension. 이제 행렬의 크기가 2000이면이 코드를 실행하는 데 147 초가 걸리는 반면, 행렬의 크기가 2048이면 447 초가 걸립니다. 그래서 차이는 없습니다. 곱셈의 수는 (2048 * 2048 * 2048) / (2000 * 2000 * 2000) = 1.073이고 타이밍의 차이는 447/147 = 3입니다. 누군가 왜 이런 일이 발생하는지 설명 할 수 있습니까? 나는 그것이 일어나지 않는 선형 적으로 확장 될 것으로 예상했다. 나는 가장 빠른 행렬 곱셈 코드를 만들려는 것이 아니라 단순히 왜 발생하는지 이해하려고 노력하고 있습니다.

사양 : AMD Opteron 듀얼 코어 노드 (2.2GHz), 2G RAM, gcc v 4.5.0

다음과 같이 컴파일 된 프로그램 gcc -O3 simple.c

나는 이것을 Intel의 icc 컴파일러에서도 실행했으며 비슷한 결과를 보았습니다.

편집하다:

댓글 / 답변에서 제안했듯이 dimension = 2060으로 코드를 실행했으며 145 초가 걸립니다.

전체 프로그램은 다음과 같습니다.

#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>

/* change dimension size as needed */
const int dimension = 2048;
struct timeval tv; 

double timestamp()
{
        double t;
        gettimeofday(&tv, NULL);
        t = tv.tv_sec + (tv.tv_usec/1000000.0);
        return t;
}

int main(int argc, char *argv[])
{
        int i, j, k;
        double *A, *B, *C, start, end;

        A = (double*)malloc(dimension*dimension*sizeof(double));
        B = (double*)malloc(dimension*dimension*sizeof(double));
        C = (double*)malloc(dimension*dimension*sizeof(double));

        srand(292);

        for(i = 0; i < dimension; i++)
                for(j = 0; j < dimension; j++)
                {   
                        A[dimension*i+j] = (rand()/(RAND_MAX + 1.0));
                        B[dimension*i+j] = (rand()/(RAND_MAX + 1.0));
                        C[dimension*i+j] = 0.0;
                }   

        start = timestamp();
        for(i = 0; i < dimension; i++)
                for(j = 0; j < dimension; j++)
                        for(k = 0; k < dimension; k++)
                                C[dimension*i+j] += A[dimension*i+k] *
                                        B[dimension*k+j];

        end = timestamp();
        printf("\nsecs:%f\n", end-start);

        free(A);
        free(B);
        free(C);

        return 0;
}

---

내 추측은 다음과 같다. 캐시

2000 double초의 2 개 행을 캐시에 넣을 수 있습니다 . 32kb L1 캐시보다 약간 적습니다. (방을 떠나면서 다른 필요한 것)

하지만 2048까지 올리면 전체 캐시를 사용합니다 (다른 것을위한 공간이 필요하기 때문에 일부 유출).

캐시 정책이 LRU라고 가정하면 캐시를 조금만 넘치면 전체 행이 반복적으로 플러시되고 L1 캐시로 다시로드됩니다.

다른 가능성은 2의 거듭 제곱으로 인한 캐시 연관성입니다. 프로세서가 2-way L1 결합이라고 생각하지만이 경우에는 중요하지 않다고 생각합니다. (하지만 어쨌든 아이디어를 던질 것입니다)

가능한 설명 2 : L2 캐시의 수퍼 정렬로 인해 충돌 캐시 누락.

귀하의 B배열은 열 반복되고있다. 따라서 접근이 제한됩니다. 총 데이터 크기는 2k x 2k매트릭스 당 약 32MB입니다. L2 캐시보다 훨씬 큽니다.

데이터가 완벽하게 정렬되지 않으면 B에 적절한 공간적 지역성을 갖게됩니다. 행을 호핑하고 캐시 라인 당 하나의 요소 만 사용하더라도 캐시 라인은 L2 캐시에 남아 중간 루프의 다음 반복에서 재사용됩니다.

그러나 데이터가 완벽하게 정렬되면 (2048 년) 이러한 홉은 모두 동일한 "캐시 방식"에 도달하고 L2 캐시 연관성을 훨씬 초과합니다. 따라서 액세스 된 캐시 라인은 B다음 반복을 위해 캐시에 남아 있지 않습니다. 대신, 그들은 숫양에서 끝까지 당겨야합니다.

---

당신은 확실히 내가 캐시 공명 이라고 부르는 것을 얻고 있습니다. 이것은 aliasing 과 유사 하지만 정확히 동일하지는 않습니다. 설명하겠습니다.

캐시는 소프트웨어의 배열과 달리 주소의 한 부분을 추출하여 테이블의 인덱스로 사용하는 하드웨어 데이터 구조입니다. (사실, 우리는 그것들을 하드웨어에서 배열이라고 부릅니다.) 캐시 배열에는 데이터의 캐시 라인과 태그가 포함되어 있습니다. 때로는 배열의 인덱스 당 하나의 항목 (직접 매핑), 때로는 여러 항목 (N-way 집합 연관성)이 있습니다. 주소의 두 번째 부분이 추출되어 배열에 저장된 태그와 비교됩니다. 인덱스와 태그는 함께 캐시 라인 메모리 주소를 고유하게 식별합니다. 마지막으로 나머지 주소 비트는 액세스 크기와 함께 캐시 라인에서 주소가 지정된 바이트를 식별합니다.

일반적으로 인덱스와 태그는 단순한 비트 필드입니다. 따라서 메모리 주소는 다음과 같습니다.

  ...Tag... | ...Index... | Offset_within_Cache_Line

(때때로 인덱스와 태그는 해시입니다. 예를 들어, 인덱스 인 중간 범위 비트에 대한 다른 비트의 XOR 몇 개가 있습니다. 훨씬 더 드물게, 때로는 인덱스 및 드물게 태그가 캐시 라인 주소 모듈로를 취하는 것과 같은 것입니다. 이보다 복잡한 인덱스 계산은 여기서 설명하는 공명 문제를 해결하려는 시도입니다. 모두가 어떤 형태의 공명을 겪지 만 가장 단순한 비트 필드 추출 방식은 일반적인 액세스 패턴에 대한 공명을 겪습니다.)

그래서, 전형적인 값들 ... "Opteron Dual Core"의 많은 다른 모델들이 있습니다. 그리고 나는 당신이 가지고있는 것을 지정하는 어떤 것도 여기에서 볼 수 없습니다. 무작위로 하나를 선택하면 AMD 웹 사이트에서 가장 최근에 볼 수있는 설명서, AMD 제품군 15h 모델 00h-0Fh 용 BKDG (Bios and Kernel Developer 's Guide , 2012 년 3 월 12 일).

(Family 15h = 가장 최근의 하이 엔드 프로세서 인 Bulldozer 제품군-BKDG는 사용자가 정확히 설명하는 제품 번호를 모르지만 듀얼 코어를 언급합니다. 그러나 어쨌든 동일한 공명 개념이 모든 프로세서에 적용됩니다. 캐시 크기 및 연관성과 같은 매개 변수가 약간 다를 수 있습니다.)

p.33에서 :

AMD 제품군 15h 프로세서에는 2 개의 128 비트 포트가있는 16KB, 4 방향 예측 L1 데이터 캐시가 포함되어 있습니다. 이것은주기 당 최대 2 개의 128 바이트로드를 지원하는 연속 기입 캐시입니다. 16 개의 뱅크로 나뉘며 각각 16 바이트 너비입니다. [...] L1 캐시의 주어진 뱅크에서 단일 주기로 하나의로드 만 수행 할 수 있습니다.

요약하자면 :

• 64 바이트 캐시 라인 => 캐시 라인 내 6 오프셋 비트

• 16KB / 4-way => 공명은 4KB입니다.

  즉, 주소 비트 0-5는 캐시 라인 오프셋입니다.

• 16KB / 64B cache lines => 2^14/2^6 = 2^8=256 cache lines in the cache.
  (Bugfix: I originally miscalculated this as 128. that I have fixed all dependencies.)

• 4 way associative => 256/4 = 64 indexes in the cache array. I (Intel) call these "sets".

  i.e. you can consider the cache to be an array of 32 entries or sets, each entry containing 4 cache lines ad their tags. (It's more complicated than this, but that's okay).

(By the way, the terms "set" and "way" have varying definitions.)

• there are 6 index bits, bits 6-11 in the simplest scheme.

  This means that any cache lines that have exactly the same values in the index bits, bits 6-11, will map to the same set of the cache.

Now look at your program.

C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];

Loop k is the innermost loop. The base type is double, 8 bytes. If dimension=2048, i.e. 2K, then successive elements of B[dimension*k+j] accessed by the loop will be 2048 * 8 = 16K bytes apart. They will all map to the same set of the L1 cache - they will all have the same index in the cache. Which means that, instead of there being 256 cache lines in the cache available for use there will only be 4 - the "4-way associativity" of the cache.

I.e. you will probably get a cache miss every 4 iterations around this loop. Not good.

(Actually, things are a little more complicated. But the above is a good first understanding. The addresses of entries of B mentioned above is a virtual address. So there might be slightly different physical addresses. Moreover, Bulldozer has a way predictive cache, probably using virtual addresses bits so that it doesn't have to wait for a virtual to physical address translation. But, in any case: your code has a "resonance" of 16K. The L1 data cache has a resonance of 16K. Not good.)]

If you change the dimension just a little bit, e.g. to 2048+1, then the addresses of array B will be spread across all of the sets of the cache. And you will get significantly fewer cache misses.

It is a fairly common optimization to pad your arrays, e.g. to change 2048 to 2049, to avoid this srt of resonance. But "cache blocking is an even more important optimization. http://suif.stanford.edu/papers/lam-asplos91.pdf

---

In addition to the cache line resonance, there are other things going on here. For example, the L1 cache has 16 banks, each 16 bytes wide. With dimension = 2048, successive B accesses in the inner loop will always go to the same bank. So they can't go in parallel - and if the A access happens to go to the same bank, you will lose.

I don't think, looking at it, that this is as big as the cache resonance.

And, yes, possibly, there may be aliasing going. E.g. the STLF (Store To Load Forwarding buffers) may be comparing only using a small bitfield, and getting false matches.

(Actually, if you think about it, resonance in the cache is like aliasing, related to the use of bitfields. Resonance is caused by multiple cache lines mapping the same set, not being spread arond. Alisaing is caused by matching based on incomplete address bits.)

---

Overall, my recommendation for tuning:

1. Try cache blocking without any further analysis. I say this because cache blocking is easy, and it is very likely that this is all you would need to do.

2. After that, use VTune or OProf. Or Cachegrind. Or ...

3. Better yet, use a well tuned library routine to do matrix multiply.

---

There are several possible explanations. One likely explanation is what Mysticial suggests: exhaustion of a limited resource (either cache or TLB). Another likely possibility is a false aliasing stall, which can occur when consecutive memory accesses are separated by a multiple of some power-of-two (often 4KB).

You can start to narrow down what's at work by plotting time/dimension^3 for a range of values. If you have blown a cache or exhausted TLB reach, you will see a more-or-less flat section followed by a sharp rise between 2000 and 2048, followed by another flat section. If you are seeing aliasing-related stalls, you will see a more-or-less flat graph with a narrow spike upward at 2048.

Of course, this has diagnostic power, but it is not conclusive. If you want to conclusively know what the source of the slowdown is, you will want to learn about performance counters, which can answer this sort of question definitively.

---

A couple of answers mentioned L2 Cache problems.

You can actually verify this with a cache simulation. Valgrind's cachegrind tool can do that.

valgrind --tool=cachegrind --cache-sim=yes your_executable

Set the command line parameters so they match with your CPU's L2 parameters.

Test it with different matrix sizes, you'll probably see a sudden increase in L2 miss ratio.

---

I know this is waaaay too old, but I'll take a bite. It's (as it has been said) a cache issue what causes the slowdown at around powers of two. But there is another problem with this: it's too slow. If you look at your compute loop.

for(i = 0; i < dimension; i++)
    for(j = 0; j < dimension; j++)
        for(k = 0; k < dimension; k++)
            C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];

The inner-most loop changes k by 1 each iteration, meaning you access just 1 double away from the last element you used of A but a whole 'dimension' doubles away from the last element of B. This doesn't take any advantage of the caching of the elements of B.

If you change this to:

for(i = 0; i < dimension; i++)
    for(j = 0; j < dimension; j++)
        for(k = 0; k < dimension; k++)
            C[dimension*i+k] += A[dimension*i+j] * B[dimension*j+k];

You get the exact same results (modulo double addition associativity errors), but it's a whole lot more cache-friendly (local). I tried it and it gives substantial improvements. This can be summarized as

Don't multiply matrices by definition, but rather, by rows

---

Example of the speed-up (I changed your code to take the dimension as an argument)

$ diff a.c b.c
42c42
<               C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];
---
>               C[dimension*i+k] += A[dimension*i+j] * B[dimension*j+k];
$ make a
cc     a.c   -o a
$ make b
cc     b.c   -o b
$ ./a 1024

secs:88.732918
$ ./b 1024

secs:12.116630

---

As a bonus (and what makes this related to this question) is that this loop doesn't suffer from the previous problem.

If you already knew all of this, then I apologize!

참고URL : https://stackoverflow.com/questions/7905760/matrix-multiplication-small-difference-in-matrix-size-large-difference-in-timi