Dlaczego transponowanie macierzy 512x512 jest znacznie wolniejsze niż transponowanie macierzy 513x513?

Po przeprowadzeniu niektórych eksperymentów na matrycach kwadratowych o różnych rozmiarach pojawił się wzór. Niezmiennie transpozycja macierzy rozmiaru 2^njest wolniejsza niż transpozycja macierzy rozmiaru2^n+1 . W przypadku małych wartości nróżnica nie jest duża.

Duże różnice występują jednak w przypadku wartości 512. (przynajmniej dla mnie)

Oświadczenie: Wiem, że funkcja nie transponuje macierzy z powodu podwójnej zamiany elementów, ale to nie robi różnicy.

Podąża za kodem:

#define SAMPLES 1000
#define MATSIZE 512

#include <time.h>
#include <iostream>
int mat[MATSIZE][MATSIZE];

void transpose()
{
   for ( int i = 0 ; i < MATSIZE ; i++ )
   for ( int j = 0 ; j < MATSIZE ; j++ )
   {
       int aux = mat[i][j];
       mat[i][j] = mat[j][i];
       mat[j][i] = aux;
   }
}

int main()
{
   //initialize matrix
   for ( int i = 0 ; i < MATSIZE ; i++ )
   for ( int j = 0 ; j < MATSIZE ; j++ )
       mat[i][j] = i+j;

   int t = clock();
   for ( int i = 0 ; i < SAMPLES ; i++ )
       transpose();
   int elapsed = clock() - t;

   std::cout << "Average for a matrix of " << MATSIZE << ": " << elapsed / SAMPLES;
}

Zmiana MATSIZEpozwala nam zmienić rozmiar (duh!). Zamieściłem dwie wersje na ideone:

• rozmiar 512 - średnio 2,46 ms - http://ideone.com/1PV7m
• rozmiar 513 - średnio 0,75 ms - http://ideone.com/NShpo

W moim środowisku (MSVS 2010, pełne optymalizacje) różnica jest podobna:

• rozmiar 512 - średnio 2,19 ms
• rozmiar 513 - średnio 0,57 ms

Dlaczego to się dzieje?

Wyjaśnienie pochodzi od Agner Fog in Optimization software in C ++ i ogranicza się do sposobu dostępu do danych i przechowywania ich w pamięci podręcznej.

Aby uzyskać warunki i szczegółowe informacje, zobacz wpis wiki dotyczący buforowania , zawęzię go tutaj.

Pamięć podręczna jest zorganizowana w zestawy i linie . Jednocześnie używany jest tylko jeden zestaw, z którego można użyć dowolnego z zawartych w nim wierszy. Pamięć, jaką może dublować linia razy liczba linii, daje nam rozmiar pamięci podręcznej.

Dla określonego adresu pamięci możemy obliczyć, który zestaw powinien go odzwierciedlać za pomocą wzoru:

set = ( address / lineSize ) % numberOfsets

Tego rodzaju formuła idealnie zapewnia równomierny rozkład między zestawami, ponieważ każdy adres pamięci jest równie prawdopodobne do odczytania ( idealnie powiedziałem ).

Oczywiste jest, że mogą się nakładać. W przypadku braku pamięci podręcznej pamięć jest odczytywana w pamięci podręcznej, a stara wartość jest zastępowana. Pamiętaj, że każdy zestaw ma pewną liczbę wierszy, z których najmniej ostatnio używany jest nadpisywany nowo odczytaną pamięcią.

Spróbuję nieco podążać za przykładem Agnera:

Załóżmy, że każdy zestaw ma 4 wiersze, każdy zawierający 64 bajty. Najpierw próbujemy odczytać adres 0x2710, który jest w zestawie 28. I wtedy też próbować odczytać adresy 0x2F00, 0x3700, 0x3F00i 0x4700. Wszystkie należą do tego samego zestawu. Przed przeczytaniem 0x4700wszystkie linie w zestawie byłyby zajęte. Odczyt tej pamięci eksmituje istniejącą linię w zestawie, linię, która początkowo trzymała 0x2710. Problem polega na tym, że czytamy adresy, które są (na przykład) 0x800osobno. To jest krytyczny krok (ponownie, w tym przykładzie).

Krytyczny krok można również obliczyć:

criticalStride = numberOfSets * lineSize

Zmienne w odstępach criticalStridelub wielokrotne odstępy rywalizują o te same linie pamięci podręcznej.

To część teoretyczna. Następnie wyjaśnienie (również Agner, uważnie je śledzę, aby uniknąć błędów):

Załóżmy, że macierz 64x64 (pamiętaj, że efekty różnią się w zależności od pamięci podręcznej) z pamięcią podręczną 8 KB, 4 liniami na zestaw * rozmiar linii 64 bajty. Każda linia może pomieścić 8 elementów w matrycy (64-bit int).

Krytyczny krok wynosiłby 2048 bajtów, co odpowiada 4 wierszom macierzy (która ma ciągłą pamięć).

Załóżmy, że przetwarzamy wiersz 28. Próbujemy pobrać elementy tego wiersza i zamienić je elementami z kolumny 28. Pierwsze 8 elementów tego wiersza tworzy wiersz pamięci podręcznej, ale przejdzie do 8 różnych buforuj linie w kolumnie 28. Pamiętaj, że krytyczny krok to 4 rzędy od siebie (4 kolejne elementy w kolumnie).

Gdy element 16 zostanie osiągnięty w kolumnie (4 wiersze pamięci podręcznej na zestaw i 4 rzędy od siebie = problem), element ex-0 zostanie usunięty z pamięci podręcznej. Gdy dotrzemy do końca kolumny, wszystkie poprzednie wiersze pamięci podręcznej zostałyby utracone i wymagałyby ponownego załadowania przy dostępie do następnego elementu (cała linia jest zastępowana).

Posiadanie rozmiaru, który nie jest wielokrotnością krytycznego kroku, zaburza ten idealny scenariusz na wypadek katastrofy, ponieważ nie mamy już do czynienia z elementami, które są krytycznie oddalone od siebie w pionie, więc liczba przeładowań pamięci podręcznej jest znacznie zmniejszona.

Kolejne zrzeczenie się odpowiedzialności - właśnie wyjaśniłem wyjaśnienie i mam nadzieję, że je przybiłem, ale mogę się mylić. W każdym razie czekam na odpowiedź (lub potwierdzenie) od Mysticial . :)

Luchian wyjaśnia, dlaczego tak się dzieje, ale pomyślałem, że dobrym pomysłem byłoby przedstawienie jednego możliwego rozwiązania tego problemu, a jednocześnie przedstawienie trochę nieświadomych algorytmów pamięci podręcznej.

Twój algorytm zasadniczo:

for (int i = 0; i < N; i++) 
   for (int j = 0; j < N; j++) 
        A[j][i] = A[i][j];

co jest po prostu okropne dla współczesnego procesora. Jednym z rozwiązań jest poznanie szczegółowych informacji o systemie pamięci podręcznej i dostosowanie algorytmu, aby uniknąć tych problemów. Działa świetnie, o ile znasz te szczegóły. Nie jest to szczególnie przenośne.

Czy możemy to zrobić lepiej? Tak, możemy: Ogólnym podejściem do tego problemu są nieświadome algorytmy pamięci podręcznej, które, jak sama nazwa wskazuje, pozwalają uniknąć zależności od konkretnych rozmiarów pamięci podręcznej [1]

Rozwiązanie wyglądałoby tak:

void recursiveTranspose(int i0, int i1, int j0, int j1) {
    int di = i1 - i0, dj = j1 - j0;
    const int LEAFSIZE = 32; // well ok caching still affects this one here
    if (di >= dj && di > LEAFSIZE) {
        int im = (i0 + i1) / 2;
        recursiveTranspose(i0, im, j0, j1);
        recursiveTranspose(im, i1, j0, j1);
    } else if (dj > LEAFSIZE) {
        int jm = (j0 + j1) / 2;
        recursiveTranspose(i0, i1, j0, jm);
        recursiveTranspose(i0, i1, jm, j1);
    } else {
    for (int i = i0; i < i1; i++ )
        for (int j = j0; j < j1; j++ )
            mat[j][i] = mat[i][j];
    }
}

Nieco bardziej złożony, ale krótki test pokazuje coś całkiem interesującego na moim starożytnym e8400 z wersją VS2010 x64, kod testowy dla MATSIZE 8192

int main() {
    LARGE_INTEGER start, end, freq;
    QueryPerformanceFrequency(&freq);
    QueryPerformanceCounter(&start);
    recursiveTranspose(0, MATSIZE, 0, MATSIZE);
    QueryPerformanceCounter(&end);
    printf("recursive: %.2fms\n", (end.QuadPart - start.QuadPart) / (double(freq.QuadPart) / 1000));

    QueryPerformanceCounter(&start);
    transpose();
    QueryPerformanceCounter(&end);
    printf("iterative: %.2fms\n", (end.QuadPart - start.QuadPart) / (double(freq.QuadPart) / 1000));
    return 0;
}

results: 
recursive: 480.58ms
iterative: 3678.46ms

Edycja: O wpływie rozmiaru: Jest znacznie mniej wyraźny, choć nadal do pewnego stopnia zauważalny, ponieważ używamy iteracyjnego rozwiązania jako węzła liścia zamiast rekurencji do 1 (zwykła optymalizacja algorytmów rekurencyjnych). Jeśli ustawimy LEAFSIZE = 1, pamięć podręczna nie ma na mnie wpływu [ 8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms- to margines błędu, wahania są w obszarze 100 ms; ten „test” nie jest czymś, z czym czułbym się zbyt dobrze, gdybyśmy chcieli całkowicie dokładnych wartości])

[1] Źródła dla tych rzeczy: Cóż, jeśli nie możesz dostać wykładu od kogoś, kto pracował z Leisersonem i współpracować w tej sprawie. Zakładam, że ich artykuły są dobrym punktem wyjścia. Algorytmy te są wciąż dość rzadko opisywane - CLR ma jeden przypis na ich temat. To wciąż świetny sposób na zaskoczenie ludzi.

Edycja (uwaga: nie jestem tym, który opublikował tę odpowiedź; chciałem ją tylko dodać):
Oto pełna wersja powyższego kodu w C ++:

template<class InIt, class OutIt>
void transpose(InIt const input, OutIt const output,
    size_t const rows, size_t const columns,
    size_t const r1 = 0, size_t const c1 = 0,
    size_t r2 = ~(size_t) 0, size_t c2 = ~(size_t) 0,
    size_t const leaf = 0x20)
{
    if (!~c2) { c2 = columns - c1; }
    if (!~r2) { r2 = rows - r1; }
    size_t const di = r2 - r1, dj = c2 - c1;
    if (di >= dj && di > leaf)
    {
        transpose(input, output, rows, columns, r1, c1, (r1 + r2) / 2, c2);
        transpose(input, output, rows, columns, (r1 + r2) / 2, c1, r2, c2);
    }
    else if (dj > leaf)
    {
        transpose(input, output, rows, columns, r1, c1, r2, (c1 + c2) / 2);
        transpose(input, output, rows, columns, r1, (c1 + c2) / 2, r2, c2);
    }
    else
    {
        for (ptrdiff_t i1 = (ptrdiff_t) r1, i2 = (ptrdiff_t) (i1 * columns);
            i1 < (ptrdiff_t) r2; ++i1, i2 += (ptrdiff_t) columns)
        {
            for (ptrdiff_t j1 = (ptrdiff_t) c1, j2 = (ptrdiff_t) (j1 * rows);
                j1 < (ptrdiff_t) c2; ++j1, j2 += (ptrdiff_t) rows)
            {
                output[j2 + i1] = input[i2 + j1];
            }
        }
    }
}

Jako ilustrację wyjaśnienia w odpowiedzi Luchiana Grigore'a , oto jak wygląda obecność pamięci podręcznej macierzy dla dwóch przypadków macierzy 64x64 i 65x65 (szczegółowe informacje na temat liczb znajdują się w powyższym linku).

Kolory w poniższych animacjach oznaczają:

• - nie w pamięci podręcznej,
• - w pamięci podręcznej,
• - trafienie w pamięć podręczną,
• - po prostu czytaj z pamięci RAM,
• - chybienie w pamięci podręcznej.

Obudowa 64x64:

Zauważ, że prawie każdy dostęp do nowego wiersza powoduje brak pamięci podręcznej. A teraz, jak wygląda normalny przypadek, macierz 65x65:

Tutaj widać, że większość dostępów po początkowej rozgrzewce to trafienia do pamięci podręcznej. W ten sposób pamięć podręczna procesora ma działać w ogóle.

_{Kod generujący ramki dla powyższych animacji można zobaczyć tutaj .}