Tengo una matriz (relativamente grande) que necesito transponer. Por ejemplo supongamos que mi matriz es
a b c d e f
g h i j k l
m n o p q r
Quiero que el resultado sea el siguiente
a g m
b h n
c I o
d j p
e k q
f l r
¿Cuál es la forma más rápida de hacerlo?
Es una buena pregunta. Hay muchas razones para querer transponer la matriz en memoria en lugar de simplemente intercambiar coordenadas, por ejemplo, en la multiplicación de matrices y en el barrido gaussiano.
En primer lugar permítanme enumerar una de las funciones que utilizo para la transposición (EDIT: por favor, consulte el final de mi respuesta donde he encontrado una solución mucho más rápida)
void transpose(float *src, float *dst, const int N, const int M) {
#pragma omp parallel for
for(int n = 0; n<N*M; n++) {
int i = n/N;
int j = n%N;
dst[n] = src[M*j + i];
}
}
Ahora vamos a ver por qué la transposición es útil. Consideremos la multiplicación de matrices C = A*B. Podríamos hacerlo de esta manera.
for(int i=0; i<N; i++) {
for(int j=0; j<K; j++) {
float tmp = 0;
for(int l=0; l<M; l++) {
tmp += A[M*i+l]*B[K*l+j];
}
C[K*i + j] = tmp;
}
}
De esa manera, sin embargo, va a tener un montón de pérdidas de caché. Una solución mucho más rápida es tomar la transposición de B primero
transpose(B);
for(int i=0; i<N; i++) {
for(int j=0; j<K; j++) {
float tmp = 0;
for(int l=0; l<M; l++) {
tmp += A[M*i+l]*B[K*j+l];
}
C[K*i + j] = tmp;
}
}
transpose(B);
La multiplicación de matrices es O(n^3) y la transposición es O(n^2), por lo que tomar la transposición debería tener un efecto insignificante en el tiempo de cálculo (para grandes n
). En la multiplicación de matrices loop tiling es aún más eficaz que tomar la transposición, pero eso es mucho más complicado.
Me gustaría saber una manera más rápida de hacer la transposición (Edición: He encontrado una solución más rápida, ver el final de mi respuesta). Cuando Haswell/AVX2 salga en unas semanas tendrá una función gather. No sé si eso será útil en este caso, pero podría imagen de reunir una columna y escribir una fila. Tal vez hará que la transposición innecesario.
Para Gaussian manchas lo que haces es manchar horizontalmente y luego manchar verticalmente. Pero manchar verticalmente tiene el problema de caché así que lo que haces es
Smear image horizontally
transpose output
Smear output horizontally
transpose output
Aquí hay un documento de Intel que explica que http://software.intel.com/en-us/articles/iir-gaussian-blur-filter-implementation-using-intel-advanced-vector-extensions
Por último, lo que hago en realidad en la multiplicación de matrices (y en el emborronamiento gaussiano) no es tomar exactamente la transposición, sino tomar la transposición en anchos de un cierto tamaño de vector (por ejemplo, 4 u 8 para SSE/AVX). Aquí está la función que utilizo
void reorder_matrix(const float* A, float* B, const int N, const int M, const int vec_size) {
#pragma omp parallel for
for(int n=0; n<M*N; n++) {
int k = vec_size*(n/N/vec_size);
int i = (n/vec_size)%N;
int j = n%vec_size;
B[n] = A[M*i + k + j];
}
}
EDIT:
He probado varias funciones para encontrar la transposición más rápida para matrices grandes. Al final, el resultado más rápido es utilizar el bloqueo de bucle con block_size=16
(Edit: He encontrado una solución más rápida utilizando SSE y el bloqueo de bucle - ver más abajo). Este código funciona para cualquier matriz NxM (es decir, la matriz no tiene que ser cuadrada).
inline void transpose_scalar_block(float *A, float *B, const int lda, const int ldb, const int block_size) {
#pragma omp parallel for
for(int i=0; i<block_size; i++) {
for(int j=0; j<block_size; j++) {
B[j*ldb + i] = A[i*lda +j];
}
}
}
inline void transpose_block(float *A, float *B, const int n, const int m, const int lda, const int ldb, const int block_size) {
#pragma omp parallel for
for(int i=0; i<n; i+=block_size) {
for(int j=0; j<m; j+=block_size) {
transpose_scalar_block(&A[i*lda +j], &B[j*ldb + i], lda, ldb, block_size);
}
}
}
Los valores lda
y ldb
son el ancho de la matriz. Tienen que ser múltiplos del tamaño del bloque. Para encontrar los valores y asignar la memoria para, por ejemplo, una matriz de 3000x1001, hago algo como esto
#define ROUND_UP(x, s) (((x)+((s)-1)) & -(s))
const int n = 3000;
const int m = 1001;
int lda = ROUND_UP(m, 16);
int ldb = ROUND_UP(n, 16);
float *A = (float*)_mm_malloc(sizeof(float)*lda*ldb, 64);
float *B = (float*)_mm_malloc(sizeof(float)*lda*ldb, 64);
Para 3000x1001 esto devuelve ldb = 3008
y lda = 1008
Edit:
He encontrado una solución aún más rápida usando intrínsecos SSE:
inline void transpose4x4_SSE(float *A, float *B, const int lda, const int ldb) {
__m128 row1 = _mm_load_ps(&A[0*lda]);
__m128 row2 = _mm_load_ps(&A[1*lda]);
__m128 row3 = _mm_load_ps(&A[2*lda]);
__m128 row4 = _mm_load_ps(&A[3*lda]);
_MM_TRANSPOSE4_PS(row1, row2, row3, row4);
_mm_store_ps(&B[0*ldb], row1);
_mm_store_ps(&B[1*ldb], row2);
_mm_store_ps(&B[2*ldb], row3);
_mm_store_ps(&B[3*ldb], row4);
}
inline void transpose_block_SSE4x4(float *A, float *B, const int n, const int m, const int lda, const int ldb ,const int block_size) {
#pragma omp parallel for
for(int i=0; i<n; i+=block_size) {
for(int j=0; j<m; j+=block_size) {
int max_i2 = i+block_size < n ? i + block_size : n;
int max_j2 = j+block_size < m ? j + block_size : m;
for(int i2=i; i2<max_i2; i2+=4) {
for(int j2=j; j2<max_j2; j2+=4) {
transpose4x4_SSE(&A[i2*lda +j2], &B[j2*ldb + i2], lda, ldb);
}
}
}
}
}
template <class T>
void transpose( std::vector< std::vector<T> > a,
std::vector< std::vector<T> > b,
int width, int height)
{
for (int i = 0; i < width; i++)
{
for (int j = 0; j < height; j++)
{
b[j][i] = a[i][j];
}
}
}