Por que a transposição de uma matriz de 512x512 é muito mais lenta que a transposição de uma matriz de 513x513?

Depois de realizar alguns experimentos em matrizes quadradas de tamanhos diferentes, surgiu um padrão. Invariavelmente, a transposição de uma matriz de tamanho 2^né mais lenta que a transposição de uma de tamanho2^n+1 . Para valores pequenos de n, a diferença não é grande.

Grandes diferenças ocorrem, no entanto, acima de um valor de 512. (pelo menos para mim)

Isenção de responsabilidade: Eu sei que a função não transpõe a matriz por causa da troca dupla de elementos, mas não faz diferença.

Segue o código:

#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;
}

Mudar MATSIZEnos permite alterar o tamanho (duh!). Postei duas versões no ideone:

• tamanho 512 - média 2,46 ms - http://ideone.com/1PV7m
• tamanho 513 - média 0,75 ms - http://ideone.com/NShpo

No meu ambiente (MSVS 2010, otimizações completas), a diferença é semelhante:

• tamanho 512 - média 2.19 ms
• tamanho 513 - média 0,57 ms

Por que isso está acontecendo?

A explicação vem de Agner Fog em software Optimizing em C ++ e reduz a forma como os dados são acessados ​​e armazenados no cache.

Para termos e informações detalhadas, consulte o entrada wiki sobre armazenamento em cache , vou reduzi-la aqui.

Um cache é organizado em conjuntos e linhas . Por vez, apenas um conjunto é usado, do qual qualquer uma das linhas que ele contém pode ser usada. A memória que uma linha pode espelhar vezes o número de linhas nos fornece o tamanho do cache.

Para um endereço de memória específico, podemos calcular qual conjunto deve espelhá-lo com a fórmula:

set = ( address / lineSize ) % numberOfsets

Idealmente, esse tipo de fórmula fornece uma distribuição uniforme entre os conjuntos, porque cada endereço de memória tem a probabilidade de ser lido (eu disse idealmente ).

É claro que podem ocorrer sobreposições. Em caso de falta de cache, a memória é lida no cache e o valor antigo é substituído. Lembre-se de que cada conjunto possui um número de linhas, das quais a menos usada recentemente é substituída pela memória recém-lida.

Vou tentar seguir um pouco o exemplo de Agner:

Suponha que cada conjunto tenha 4 linhas, cada uma contendo 64 bytes. Primeiro tentamos ler o endereço 0x2710, que entra em conjunto 28. E então nós também tentar ler endereços 0x2F00, 0x3700, 0x3F00e 0x4700. Todos esses pertencem ao mesmo conjunto. Antes da leitura 0x4700, todas as linhas do conjunto teriam sido ocupadas. A leitura dessa memória elimina uma linha existente no conjunto, a linha que inicialmente estava mantendo 0x2710. O problema está no fato de lermos endereços que são (neste exemplo) 0x800separados. Este é o passo crítico (novamente, para este exemplo).

O passo crítico também pode ser calculado:

criticalStride = numberOfSets * lineSize

Variáveis ​​espaçadas criticalStrideou um múltiplo separado disputam as mesmas linhas de cache.

Esta é a parte da teoria. A seguir, a explicação (também Agner, estou acompanhando de perto para evitar erros):

Suponha uma matriz de 64x64 (lembre-se, os efeitos variam de acordo com o cache) com um cache de 8kb, 4 linhas por conjunto * tamanho de linha de 64 bytes. Cada linha pode conter 8 dos elementos na matriz (64 bits int).

O passo crítico seria 2048 bytes, que correspondem a 4 linhas da matriz (que é contínua na memória).

Suponha que estamos processando a linha 28. Estamos tentando pegar os elementos desta linha e trocá-los pelos elementos da coluna 28. Os primeiros 8 elementos da linha formam uma linha de cache, mas eles entram em 8 diferentes linhas de cache na coluna 28. Lembre-se de que o passo crítico está separado por 4 linhas (4 elementos consecutivos em uma coluna).

Quando o elemento 16 for alcançado na coluna (4 linhas de cache por conjunto e 4 linhas separadas = problema), o elemento ex-0 será removido do cache. Quando chegamos ao final da coluna, todas as linhas de cache anteriores teriam sido perdidas e necessárias para recarregar o acesso ao próximo elemento (a linha inteira é substituída).

Ter um tamanho que não seja múltiplo do passo crítico atrapalha esse cenário perfeito para um desastre, já que não estamos mais lidando com elementos que estão separados do ponto crítico na vertical, portanto o número de recargas de cache é severamente reduzido.

Outro aviso - acabei de entender a explicação e espero ter acertado em cheio, mas posso estar enganado. De qualquer forma, estou esperando por uma resposta (ou confirmação) de Mysticial . :)

Luchian dá uma explicação de por que esse comportamento ocorre, mas achei que seria uma boa idéia mostrar uma solução possível para esse problema e, ao mesmo tempo, mostrar um pouco sobre algoritmos alheios ao cache.

Seu algoritmo basicamente faz:

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

o que é simplesmente horrível para uma CPU moderna. Uma solução é conhecer os detalhes sobre o seu sistema de cache e ajustar o algoritmo para evitar esses problemas. Funciona muito bem desde que você conheça esses detalhes. Não é especialmente portátil.

Podemos fazer melhor que isso? Sim, podemos: Uma abordagem geral para esse problema são algoritmos alheios ao cache que, como o nome indica, evitam depender de tamanhos específicos de cache [1]

A solução seria assim:

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];
    }
}

Um pouco mais complexo, mas um pequeno teste mostra algo bastante interessante no meu antigo e8400 com o lançamento do VS2010 x64, código de teste para 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

Edit: Sobre a influência do tamanho: é muito menos pronunciado, embora ainda perceptível até certo ponto, é porque estamos usando a solução iterativa como um nó folha em vez de recursar até 1 (a otimização usual para algoritmos recursivos). Se definirmos LEAFSIZE = 1, o cache não terá influência para mim [ 8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms- isso está dentro da margem de erro, as flutuações estão na área de 100ms; esse "benchmark" não é algo com o qual eu me sentiria muito confortável se quiséssemos valores completamente precisos])

[1] Fontes para essas coisas: Bem, se você não conseguir uma palestra de alguém que trabalhou com Leiserson e co-nisto ... presumo que os trabalhos deles sejam um bom ponto de partida. Esses algoritmos ainda são raramente descritos - o CLR tem uma única nota de rodapé sobre eles. Ainda é uma ótima maneira de surpreender as pessoas.

Editar (nota: não fui eu quem postou esta resposta; eu só queria adicionar esta):
Aqui está uma versão completa em C ++ do código acima:

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];
            }
        }
    }
}

Como uma ilustração da explicação na resposta de Luchian Grigore , veja como é a presença do cache de matriz nos dois casos de matrizes 64x64 e 65x65 (veja o link acima para obter detalhes sobre números).

As cores nas animações abaixo significam o seguinte:

• - não está em cache,
• - em cache,
• - acerto no cache,
• - apenas leia da RAM,
• - falta de cache.

O caso de 64x64:

Observe como quase todo acesso a uma nova linha resulta em uma falha de cache. E agora, como parece o caso normal, uma matriz de 65x65:

Aqui você pode ver que a maioria dos acessos após o aquecimento inicial são ocorrências de cache. É assim que o cache da CPU deve funcionar em geral.

_{O código que gerou quadros para as animações acima pode ser visto aqui .}