Utilización de la caché en la transposición de la matriz en c

Question

Este código transpone una matriz de cuatro maneras. El primero hace escrituras secuenciales, lecturas no secuenciales. El segundo es todo lo contrario. Los siguientes dos son iguales, pero con la caché omitiendo las escrituras. Lo que parece suceder es que las escrituras secuenciales son más rápidas, y omitir el caché es más rápido. Lo que no entiendo es que, si se está omitiendo el caché, ¿por qué las escrituras secuenciales son aún más rápidas?

QueryPerformanceCounter(&before); 
for (i = 0; i < N; ++i) 
    for (j = 0; j < N; ++j) 
     tmp[i][j] = mul2[j][i]; 
QueryPerformanceCounter(&after); 
printf("Transpose 1:\t%ld\n", after.QuadPart - before.QuadPart); 

QueryPerformanceCounter(&before); 
for (j = 0; j < N; ++j) 
    for (i = 0; i < N; ++i) 
    tmp[i][j] = mul2[j][i]; 
QueryPerformanceCounter(&after); 
printf("Transpose 2:\t%ld\n", after.QuadPart - before.QuadPart); 

QueryPerformanceCounter(&before); 
for (i = 0; i < N; ++i) 
    for (j = 0; j < N; ++j) 
     _mm_stream_si32(&tmp[i][j], mul2[j][i]); 
QueryPerformanceCounter(&after); 
printf("Transpose 3:\t%ld\n", after.QuadPart - before.QuadPart); 

QueryPerformanceCounter(&before); 
for (j = 0; j < N; ++j) 
    for (i = 0; i < N; ++i) 
     _mm_stream_si32(&tmp[i][j], mul2[j][i]); 
QueryPerformanceCounter(&after); 
printf("Transpose 4:\t%ld\n", after.QuadPart - before.QuadPart); 

EDIT: La salida es

Transpose 1: 47603 
Transpose 2: 92449 
Transpose 3: 38340 
Transpose 4: 69597 

Answer 1

CPU tiene una memoria intermedia combinación de escritura para combinar escribe en una línea de caché a suceder en una ráfaga. En este caso (caché omitida para escrituras secuenciales), este buffer de combinación de escritura actúa como como un caché de una línea que hace que los resultados sean muy similares a la caché que no se salta.

Para ser exactos, en caso de que se omita la caché, las escrituras todavía están sucediendo en ráfagas en la memoria.

Comportamiento write-combining logic aquí.

Answer 2

Puede probar el diseño de memoria no lineal para que la matriz mejore la utilización de la memoria caché. Con los mosaicos flotantes 4x4 de 32 bits uno podría hacer la transposición con solo acceso único a cada línea de caché. Además, como un mosaico de bonificación, el cambio se puede hacer fácilmente con _MM_TRANSPOSE4_PS.

La transposición de una matriz muy grande sigue siendo una operación con mucha memoria. Todavía tendrá un gran ancho de banda, pero al menos la carga de la palabra del caché es casi óptima. No sé si el rendimiento podría optimizarse. Mis pruebas muestran que una laptop de pocos años logra transportar 16k * 16k (memoria 1G) en aproximadamente 300ms.

Intenté usar también _mm_stream_sd pero en realidad empeora el rendimiento por alguna razón. No entiendo la memoria no temporal escribe lo suficiente como para tener una suposición práctica de por qué el rendimiento caería con _mm_stream_ps. La razón posible es, por supuesto, que la línea de caché ya está en caché L1 lista para la operación de escritura.

Pero la parte realmente importante con la matriz no lineal podría evitar la transposición completamente y ejecutar simplemente la multiplicación en el orden amistoso. Pero solo tengo el código de transposición que estoy usando para mejorar mi conocimiento sobre la administración de memoria caché en algoritmos.

Aún no he intentado probar si la captación previa mejoraría el uso del ancho de banda de la memoria. El código actual funciona a aproximadamente 0.5 instrucciones por ciclo (un buen código amigable de caché funciona alrededor de 2 ins por ciclo en esta CPU) que deja muchos ciclos libres para las instrucciones de precaptura permitiendo incluso cálculos bastante complejos para optimizar el tiempo de captación previa en tiempo de ejecución.

código de ejemplo de mi prueba de benchmark transpose sigue.

#define MATSIZE 16384 
#define align(val, a) (val + (a - val % a)) 
#define tilewidth 4 
typedef int matrix[align(MATSIZE,tilewidth)*MATSIZE] __attribute__((aligned(64))); 


float &index(matrix m, unsigned i, unsigned j) 
{ 
    /* tiled address calculation */ 
    /* single cache line is used for 4x4 sub matrices (64 bytes = 4*4*sizeof(int) */ 
    /* tiles are arranged linearly from top to bottom */ 
    /* 
    * eg: 16x16 matrix tile positions: 
    * t1 t5 t9 t13 
    * t2 t6 t10 t14 
    * t3 t7 t11 t15 
    * t4 t8 t12 t16 
    */ 
    const unsigned tilestride = tilewidth * MATSIZE; 
    const unsigned comp0 = i % tilewidth; /* i inside tile is least significant part */ 
    const unsigned comp1 = j * tilewidth; /* next part is j multiplied by tile width */ 
    const unsigned comp2 = i/tilewidth * tilestride; 
    const unsigned add = comp0 + comp1 + comp2; 
    return m[add]; 
} 

/* Get start of tile reference */ 
float &tile(matrix m, unsigned i, unsigned j) 
{ 
    const unsigned tilestride = tilewidth * MATSIZE; 
    const unsigned comp1 = j * tilewidth; /* next part is j multiplied by tile width */ 
    const unsigned comp2 = i/tilewidth * tilestride; 
    return m[comp1 + comp2]; 

} 
template<bool diagonal> 
static void doswap(matrix mat, unsigned i, unsigned j) 
{ 
     /* special path to swap whole tile at once */ 
     union { 
      float *fs; 
      __m128 *mm; 
     } src, dst; 
     src.fs = &tile(mat, i, j); 
     dst.fs = &tile(mat, j, i); 
     if (!diagonal) { 
      __m128 srcrow0 = src.mm[0]; 
      __m128 srcrow1 = src.mm[1]; 
      __m128 srcrow2 = src.mm[2]; 
      __m128 srcrow3 = src.mm[3]; 
      __m128 dstrow0 = dst.mm[0]; 
      __m128 dstrow1 = dst.mm[1]; 
      __m128 dstrow2 = dst.mm[2]; 
      __m128 dstrow3 = dst.mm[3]; 

      _MM_TRANSPOSE4_PS(srcrow0, srcrow1, srcrow2, srcrow3); 
      _MM_TRANSPOSE4_PS(dstrow0, dstrow1, dstrow2, dstrow3); 

#if STREAMWRITE == 1 
      _mm_stream_ps(src.fs + 0, dstrow0); 
      _mm_stream_ps(src.fs + 4, dstrow1); 
      _mm_stream_ps(src.fs + 8, dstrow2); 
      _mm_stream_ps(src.fs + 12, dstrow3); 
      _mm_stream_ps(dst.fs + 0, srcrow0); 
      _mm_stream_ps(dst.fs + 4, srcrow1); 
      _mm_stream_ps(dst.fs + 8, srcrow2); 
      _mm_stream_ps(dst.fs + 12, srcrow3); 
#else 
      src.mm[0] = dstrow0; 
      src.mm[1] = dstrow1; 
      src.mm[2] = dstrow2; 
      src.mm[3] = dstrow3; 
      dst.mm[0] = srcrow0; 
      dst.mm[1] = srcrow1; 
      dst.mm[2] = srcrow2; 
      dst.mm[3] = srcrow3; 
#endif 
     } else { 
      __m128 srcrow0 = src.mm[0]; 
      __m128 srcrow1 = src.mm[1]; 
      __m128 srcrow2 = src.mm[2]; 
      __m128 srcrow3 = src.mm[3]; 

      _MM_TRANSPOSE4_PS(srcrow0, srcrow1, srcrow2, srcrow3); 

#if STREAMWRITE == 1 
      _mm_stream_ps(src.fs + 0, srcrow0); 
      _mm_stream_ps(src.fs + 4, srcrow1); 
      _mm_stream_ps(src.fs + 8, srcrow2); 
      _mm_stream_ps(src.fs + 12, srcrow3); 
#else 
      src.mm[0] = srcrow0; 
      src.mm[1] = srcrow1; 
      src.mm[2] = srcrow2; 
      src.mm[3] = srcrow3; 
#endif 
     } 
    } 
} 

static void transpose(matrix mat) 
{ 
    const unsigned xstep = 256; 
    const unsigned ystep = 256; 
    const unsigned istep = 4; 
    const unsigned jstep = 4; 
    unsigned x1, y1, i, j; 
    /* need to increment x check for y limit to allow unrolled inner loop 
    * access entries close to diagonal axis 
    */ 
    for (x1 = 0; x1 < MATSIZE - xstep + 1 && MATSIZE > xstep && xstep; x1 += xstep) 
     for (y1 = 0; y1 < std::min(MATSIZE - ystep + 1, x1 + 1); y1 += ystep) 
      for (i = x1 ; i < x1 + xstep; i += istep) { 
       for (j = y1 ; j < std::min(y1 + ystep, i); j+= jstep) 
       { 
        doswap<false>(mat, i, j); 
       } 
       if (i == j && j < (y1 + ystep)) 
        doswap<true>(mat, i, j); 
      } 

    for (i = 0 ; i < x1; i += istep) { 
     for (j = y1 ; j < std::min(MATSIZE - jstep + 1, i); j+= jstep) 
     { 
      doswap<false>(mat, i, j); 
     } 
     if (i == j) 
      doswap<true>(mat, i, j); 
    } 
    for (i = x1 ; i < MATSIZE - istep + 1; i += istep) { 
     for (j = y1 ; j < std::min(MATSIZE - jstep + 1, i); j+= jstep) 
     { 
      doswap<false>(mat, i, j); 
     } 
     if (i == j) 
      doswap<true>(mat, i, j); 
    } 
    x1 = MATSIZE - MATSIZE % istep; 
    y1 = MATSIZE - MATSIZE % jstep; 

    for (i = x1 ; i < MATSIZE; i++) 
     for (j = 0 ; j < std::min((unsigned)MATSIZE, i); j++) 
        std::swap(index(mat, i, j+0), index(mat, j+0, i)); 

    for (i = 0; i < x1; i++) 
     for (j = y1 ; j < std::min((unsigned)MATSIZE, i) ; j++) 
        std::swap(index(mat, i, j+0), index(mat, j+0, i)); 
}