Warum ist das transponieren einer matrix von 512x512 viel langsamer als transponieren einer matrix von 513x513?

Luchian gibt eine Erklärung warum dieses Verhalten passiert, aber ich dachte, es wäre eine nette Idee, um zu zeigen, eine mögliche Lösung für dieses problem und gleichzeitig zeigen Sie ein wenig über cache-oblivious-algorithmen.

Ihren Algorithmus im Grunde genommen:

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

ist einfach nur schrecklich für eine moderne CPU. Eine Lösung ist, um zu wissen, die details über Ihre cache-system und die Anpassung des Algorithmus zu vermeiden diese Probleme. Funktioniert sehr gut, solange Sie wissen, diese details.. nicht besonders portable.

Können wir besser machen? Ja, wir können: Einen Allgemeinen Ansatz für dieses problem sind cache-oblivious-algorithmen, die, wie der name schon sagt, vermeidet es, abhängig von der spezifischen cache-Größen [1]

Die Lösung würde wie folgt Aussehen:

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

Etwas komplexer, aber ein kurzer test zeigt etwas sehr Interessantes auf meinem alten e8400 mit VS2010 x64 Version, testcode für 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: Über den Einfluss der Größe: Es ist viel weniger ausgeprägt ist zwar noch spürbar um einige Grad, das ist, weil wir mit der iterativen Lösung als Blattknoten statt recursing bis zu 1 (die übliche Optimierung für rekursive algorithmen). Wenn wir LEAFSIZE = 1, wird der cache hat keinen Einfluss für mich [8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms - das ist innerhalb der Fehlergrenze, die Schwankungen sind im 100ms-Bereich, das "benchmark" ist nicht etwas, dass ich zu bequem mit, wenn wir wollten, ganz genaue Werte])

[1] Quellen für diese Sachen: Gut, wenn Sie nicht bekommen kann einen Vortrag von jemandem, der arbeitete mit Leiserson und co auf dieser.. ich nehme an, Ihre Papiere ein guter Ausgangspunkt. Diese algorithmen sind noch Recht selten beschrieben - CLR hat eine einzige Fußnote über Sie. Noch, es ist eine großartige Möglichkeit, um die Leute zu überraschen.

Bearbeiten (Hinweis: ich bin nicht derjenige, der dies geschrieben, Antwort; ich wollte nur hinzufügen das):

Hier ist eine komplette C++ - version des obigen Codes:

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