gcc -O0 outperforming -O3 on matrix sizes that are powers of 2 (matrix transpositions)

Question 1

The periodic increase of execution time must be due to the cache being only N-way associative instead of fully associative. You are witnessing hash collision related to cache line selection algorithm.

The fastest L1 cache has a smaller number of cache lines than the next level L2. In each level each cache line can be filled only from a limited set of sources.

Typical HW implementations of cache line selection algorithms will just use few bits from the memory address to determine in which cache slot the data should be written -- in HW bit shifts are free.

This causes a competition between memory ranges e.g. between addresses 0x300010 and 0x341010. In fully sequential algorithm this doesn't matter -- N is large enough for practically all algorithms of the form:

 for (i=0;i<1000;i++) a[i] += b[i] * c[i] + d[i];

But when the number of the inputs (or outputs) gets larger, which happens internally when the algorithm is optimized, having one input in the cache forces another input out of the cache.

 // one possible method of optimization with 2 outputs and 6 inputs
 // with two unrelated execution paths -- should be faster, but maybe it isn't
 for (i=0;i<500;i++) { 
       a[i]     += b[i]     * c[i]     + d[i];
       a[i+500] += b[i+500] * c[i+500] + d[i+500];
 }

A graph in Example 5: Cache Associativity illustrates 512 byte offset between matrix lines being a global worst case dimension for the particular system. When this is known, a working mitigation is to over-allocate the matrix horizontally to some other dimension char matrix[512][512 + 64].

Question 2

The improvement in performance is likely related to CPU/RAM caching.

When the data is not a power of 2, a cache line load (like 16, 32, or 64 words) transfers more than the data that is required tying up the bus—uselessly as it turns out. For a data set which is a power of 2, all of the pre-fetched data is used.

I bet if you were to disable L1 and L2 caching, the performance would be completely smooth and predictable. But it would be much slower. Caching really helps performance!

Question 3

Comment with code: In the -O3 case, with

#include <cstdlib>

extern void transpose(const size_t n, double* a)
{
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = i + 1; j < n; ++j) {
            std::swap(a[i * n + j], a[j * n + i]); // or your expanded version.
        }
    }
}

compiling with

$ g++ --version
g++ (Ubuntu/Linaro 4.8.1-10ubuntu9) 4.8.1
...
$ g++ -g1 -std=c++11 -Wall -o test.S -S test.cpp -O3

I get

_Z9transposemPd:
.LFB68:
    .cfi_startproc
.LBB2:
    testq   %rdi, %rdi
    je  .L1
    leaq    8(,%rdi,8), %r10
    xorl    %r8d, %r8d
.LBB3:
    addq    $1, %r8
    leaq    -8(%r10), %rcx
    cmpq    %rdi, %r8
    leaq    (%rsi,%rcx), %r9
    je  .L1
    .p2align 4,,10
    .p2align 3
.L10:
    movq    %r9, %rdx
    movq    %r8, %rax
    .p2align 4,,10
    .p2align 3
.L5:
.LBB4:
    movsd   (%rdx), %xmm1
    movsd   (%rsi,%rax,8), %xmm0
    movsd   %xmm1, (%rsi,%rax,8)
.LBE4:
    addq    $1, %rax
.LBB5:
    movsd   %xmm0, (%rdx)
    addq    %rcx, %rdx
.LBE5:
    cmpq    %rdi, %rax
    jne .L5
    addq    $1, %r8
    addq    %r10, %r9
    addq    %rcx, %rsi
    cmpq    %rdi, %r8
    jne .L10
.L1:
    rep ret
.LBE3:
.LBE2:
    .cfi_endproc

And something quite different if I add -m32.

(Note: it makes no difference to the assembly whether I use std::swap or your variant)

In order to understand what is causing the spikes, though, you probably want to visualize the memory operations going on.

Question 4

To add to others: g++ -std=c++11 -march=core2 -O3 -c -S - gcc version 4.8.2 (MacPorts gcc48 4.8.2_0) - x86_64-apple-darwin13.0.0 :

__Z9transposemPd:
LFB0:
        testq   %rdi, %rdi
        je      L1
        leaq    8(,%rdi,8), %r10
        xorl    %r8d, %r8d
        leaq    -8(%r10), %rcx
        addq    $1, %r8
        leaq    (%rsi,%rcx), %r9
        cmpq    %rdi, %r8
        je      L1
        .align 4,0x90
L10:
        movq    %r9, %rdx
        movq    %r8, %rax
        .align 4,0x90
L5:
        movsd   (%rdx), %xmm0
        movsd   (%rsi,%rax,8), %xmm1
        movsd   %xmm0, (%rsi,%rax,8)
        addq    $1, %rax
        movsd   %xmm1, (%rdx)
        addq    %rcx, %rdx
        cmpq    %rdi, %rax
        jne     L5
        addq    $1, %r8
        addq    %r10, %r9
        addq    %rcx, %rsi
        cmpq    %rdi, %r8
        jne     L10
L1:
        rep; ret

Basically the same as @ksfone's code, for:

#include <cstddef>

void transpose(const size_t _n, double* _A) {
    for(size_t i=0; i < _n; ++i) {
        for(size_t j=i+1; j < _n; ++j) {
            double tmp  = _A[i*_n+j];
            _A[i*_n+j] = _A[j*_n+i];
            _A[j*_n+i] = tmp;
        }
    }
}

Apart from the Mach-O 'as' differences (extra underscore, align and DWARF locations), it's the same. But very different from the OP's assembly output. A much 'tighter' inner loop.