Risultati errati durante l'uso di instrinsics SSE in Visual Studio 2010/2012 e nella modalità di rilascio

Question

Sto calcolando la media e la varianza di un array utilizzando gli intrinseci SSE. Fondamentalmente, questa è la somma dei valori e le sue piazze che può essere illustrato dal seguente programma:

int main(int argc, const char* argv[]) 
{ 
    union u 
    { 
     __m128 m; 
     float f[4]; 
    } x; 

    // Allocate memory and initialize data: [1,2,3,...stSize+1] 
    const size_t stSize = 1024; 
    float *pData = (float*) _aligned_malloc(stSize*sizeof(float), 32); 
    for (size_t s = 0; s < stSize; ++s) { 
     pData[s] = s+1; 
    } 

    // Sum and sum of squares 
    { 
     // Accumlation using SSE intrinsics 
     __m128 mEX = _mm_set_ps1(0.f); 
     __m128 mEXX = _mm_set_ps1(0.f); 
     for (size_t s = 0; s < stSize; s+=4) 
     { 
      __m128 m = _mm_load_ps(pData + s);  
      mEX = _mm_add_ps(mEX, m); 
      mEXX = _mm_add_ps(mEXX, _mm_mul_ps(m,m)); 
     } 

     // Final reduction 
     x.m = mEX; 
     double dEX = x.f[0] + x.f[1] + x.f[2] + x.f[3]; 
     x.m = mEXX; 
     double dEXX = x.f[0] + x.f[1] + x.f[2] + x.f[3]; 

     std::cout << "Sum expected: " << (stSize * stSize + stSize)/2 << std::endl; 
     std::cout << "EX: " << dEX << std::endl; 
     std::cout << "Sum of squares expected: " << 1.0/6.0 * stSize * (stSize + 1) * (2 * stSize + 1) << std::endl; 
     std::cout << "EXX: " << dEXX << std::endl; 
    } 

    // Clean up 
    _aligned_free(pData); 
} 

Ora, quando compila e si esegue il programma in modalità debug ottengo il seguente (e corretta) output:

Sum expected: 524800 
EX: 524800 
Sum of squares expected: 3.58438e+008 
EXX: 3.58438e+008 

Tuttavia, compilare ed eseguire il programma in modalità di rilascio dei seguenti (e sbagliato) i risultati vengono prodotti:

Sum expected: 524800 
EX: 524800 
Sum of squares expected: 3.58438e+008 
EXX: 3.49272e+012 

Modifica dell'ordine di accumulazione, vale a dire . EXX sia aggiornato prima EX, i risultati sono OK:

Sum expected: 524800 
EX: 524800 
Sum of squares expected: 3.58438e+008 
EXX: 3.58438e+008 

Sembra un'ottimizzazione del compilatore 'controproducente' o perché è l'ordine di esecuzione in questione? È un bug noto?

MODIFICA: Ho appena visto l'output dell'assembler. Ecco cosa ricevo (solo le parti rilevanti). Per la build di rilascio con /arch:AVX flag di compilazione abbiamo:

; 69 : // Second test: sum and sum of squares 
; 70 : { 
; 71 :  __m128 mEX = _mm_set_ps1(0.f); 
vmovaps xmm1, XMMWORD PTR __xmm@0 
mov ecx, 256    ; 00000100H 

; 72 :  __m128 mEXX = _mm_set_ps1(0.f); 
vmovaps xmm2, xmm1 
npad 12 
$LL3@main: 

; 73 :  for (size_t s = 0; s < stSize; s+=4) 
; 74 :  { 
; 75 :   __m128 m = _mm_load_ps(pData + s);  
vmovaps xmm0, xmm1 

; 76 :   mEX = _mm_add_ps(mEX, m); 
vaddps xmm1, xmm1, XMMWORD PTR [rax] 
add rax, 16 

; 77 :   mEXX = _mm_add_ps(mEXX, _mm_mul_ps(m,m)); 
vmulps xmm0, xmm0, xmm0 
vaddps xmm2, xmm0, xmm2 
dec rcx 
jne SHORT $LL3@main 

Questo è chiaramente sbagliato come questa (1) salva il risultato accumulato EX (xmm1) in xmm0 (2) accumula EX con il valore corrente (XMMWORD PTR [rax]) e (3) si accumula in EXX (xmm2) il quadrato del risultato EX accumulato precedentemente salvato in xmm0.

Al contrario, la versione senza il /arch:AVX guarda bene e come previsto:

; 69 : // Second test: sum and sum of squares 
; 70 : { 
; 71 :  __m128 mEX = _mm_set_ps1(0.f); 
movaps xmm1, XMMWORD PTR __xmm@0 
mov ecx, 256    ; 00000100H 

; 72 :  __m128 mEXX = _mm_set_ps1(0.f); 
movaps xmm2, xmm1 
npad 10 
$LL3@main: 

; 73 :  for (size_t s = 0; s < stSize; s+=4) 
; 74 :  { 
; 75 :   __m128 m = _mm_load_ps(pData + s);  
movaps xmm0, XMMWORD PTR [rax] 
add rax, 16 
dec rcx 

; 76 :   mEX = _mm_add_ps(mEX, m); 
addps xmm1, xmm0 

; 77 :   mEXX = _mm_add_ps(mEXX, _mm_mul_ps(m,m)); 
mulps xmm0, xmm0 
addps xmm2, xmm0 
jne SHORT $LL3@main 

Questo sembra davvero come un insetto. Qualcuno può confermare o confutare questo problema con una versione del compilatore diversa? (Attualmente non ho il permesso di aggiornare il compilatore)

Answer 1

Invece di eseguire manualmente l'aggiunta orizzontale, mi consiglia di utilizzare il corrispondente di istruzioni SSE _mm_hadd_ps

// Final reduction 
__m128 sum1 = _mm_hadd_ps(mEX, mEXX); 
     // == {EX[0]+EX[1], EX[2]+EX[3], EXX[0]+EXX[1], EXX[2]+EXX[3]} 
// final sum and conversion to double: 
__m128d sum2 = _mm_cvtps_pd(_mm_hadd_ps(sum1, sum1)); 
// result vector: 
double dEX_EXX[2]; // (I don't know MSVC syntax for stack aligned arrays) 
// store register to stack: (should be _mm_store_pd, if the array is aligned) 
_mm_storeu_pd(dEX_EXX, sum2); 
std::cout << "EX: " << dEX_EXX[0] << "\nEXX: " << dEX_EXX[1] << std::endl;