¿Cuál es la forma más rápida de transponer una matriz en C ++?

81

Tengo una matriz (relativamente grande) que necesito transponer. Por ejemplo, suponga 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 hacer esto?

mans
fuente
2
Eso se llama "transposición". Girar 90 grados es una noción completamente diferente.
Andy Prowl
35
Y la forma más rápida es no rotarlo, sino simplemente intercambiar el orden del índice cuando acceda a la matriz.
High Performance Mark
2
No importa qué tan rápido sea, debe acceder a todos los elementos de la matriz de todos modos.
taocp
10
@HighPerformanceMark: Supongo que depende, si luego desea acceder a la matriz repetidamente en el orden de las filas, tener una bandera "transpuesta" le afectará mucho.
Matthieu M.
3
La transposición de matrices es notoria por los problemas que causa con los cachés de memoria. Si su matriz es lo suficientemente grande como para que el rendimiento de una transposición sea significativo, y no puede evitar la transposición simplemente proporcionando una interfaz con índices intercambiados, entonces su mejor opción es usar una rutina de biblioteca existente para transponer matrices grandes. Los expertos ya han realizado este trabajo y debería utilizarlo.
Eric Postpischil

Respuestas:

131

Esta es una buena pregunta. Hay muchas razones por las que desearía realmente transponer la matriz en la memoria en lugar de simplemente intercambiar coordenadas, por ejemplo, en la multiplicación de matrices y la difusión de Gauss.

Primero permítanme enumerar una de las funciones que uso para la transposición ( EDITAR: vea el final de mi respuesta donde encontré 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 veamos por qué es útil la transposición. Considere 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 forma, sin embargo, habrá muchas pérdidas de caché. Una solución mucho más rápida es tomar primero la transposición de B

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 el ciclo de multiplicación de matrices, el mosaico es incluso más efectivo que tomar la transposición, pero eso es mucho más complicado.

Ojalá supiera una forma más rápida de hacer la transposición ( Editar: encontré una solución más rápida, vea el final de mi respuesta ). Cuando salga Haswell / AVX2 en unas pocas semanas, tendrá una función de recopilación. No sé si eso será útil en este caso, pero podría imaginarme reuniendo una columna y escribiendo una fila. Tal vez haga innecesaria la transposición.

Para el difuminado gaussiano, lo que se hace es untar horizontalmente y luego untar verticalmente. Pero manchar verticalmente tiene el problema de la caché, por lo que

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 realmente hago en la multiplicación de matrices (y en la difuminación gaussiana) 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 uso

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

EDITAR:

Probé varias funciones para encontrar la transposición más rápida para matrices grandes. Al final, el resultado más rápido es usar el bloqueo de bucle con block_size=16( Editar: encontré una solución más rápida usando SSE y 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 lday ldbson el ancho de la matriz. Estos deben 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 = 3008y lda = 1008

Editar:

Encontré 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);
                }
            }
        }
    }
}
Cole Johnson
fuente
1
Buen tiro, pero no estoy seguro de que 'La multiplicación de matrices es O (n ^ 3)', creo que es O (n ^ 2).
ulyssis2
2
@ ulyssis2 Es O (n ^ 3), a menos que use la Multiplicación de matrices de Strassen (O (n ^ 2.8074)). user2088790: Esto está muy bien hecho. Manteniendo esto en mi colección personal. :)
saurabheights
10
En caso de que alguien quiera saber quién escribió esta respuesta, fui yo. Dejé SO una vez, lo superé y regresé.
Bosón Z
1
@ ulyssis2 La multiplicación de matrices ingenua es definitivamente O (n ^ 3) y, hasta donde yo sé, los núcleos de cálculo implementan el algoritmo ingenuo (creo que esto se debe a que Strassen termina haciendo muchas más operaciones (adiciones), lo cual es malo si puedes hacer productos rápidos, pero podría estar equivocado). Es un problema abierto si la multiplicación de matrices puede ser O (n ^ 2) o no.
étale-cohomology
Por lo general, es una mejor opción confiar en una biblioteca de álgebra lineal para hacer el trabajo por usted. Las bibliotecas modernas como Intel MKL, OpenBLAS, etc.proporcionan un despacho dinámico de CPU que selecciona la mejor implementación disponible para su hardware (por ejemplo, pueden estar disponibles registros vectoriales más amplios que SSE: AVX AVX2, AVX512 ...), por lo que no No es necesario crear un programa no portátil para obtener un programa rápido.
Jorge Bellon
39

Esto dependerá de su aplicación, pero en general la forma más rápida de transponer una matriz sería invertir sus coordenadas cuando realiza una búsqueda, entonces no tiene que mover ningún dato.

Shafik Yaghmour
fuente
32
Esto es genial si es una matriz pequeña o solo la lee una vez. Sin embargo, si la matriz transpuesta es grande y necesita ser reutilizada muchas veces, aún puede guardar una versión transpuesta rápida para obtener un mejor patrón de acceso a la memoria. (+1, por cierto)
Agentlien
2
@Agentlien: ¿Por qué A [j] [i] sería más lento que A [i] [j]?
vaso de precipitados
32
@beaker Si tiene una matriz grande, diferentes filas / columnas pueden ocupar diferentes líneas / páginas de caché. En este caso, querrá iterar sobre elementos de tal manera que acceda a elementos adyacentes uno después del otro. De lo contrario, puede hacer que el acceso a cada elemento se convierta en una falta de caché, lo que destruye por completo el rendimiento.
Agentlien
10
@beaker: tiene que ver con el almacenamiento en caché a nivel de la CPU (suponiendo que la matriz es una gran cantidad de memoria), las líneas de caché son entonces líneas efectivas de la matriz y el prefetcher puede buscar las siguientes líneas. Si cambia el acceso, el caché / prefetcher de la CPU aún funciona línea por línea mientras accede columna por columna, la caída del rendimiento puede ser dramática.
Matthieu M.
2
@taocp Básicamente, necesitaría algún tipo de bandera para indicar que está transpuesta y luego la solicitud de decir (i,j)se (j,i)
mapeará
5

Algunos detalles sobre la transposición de matrices flotantes cuadradas de 4x4 (discutiré el entero de 32 bits más adelante) con hardware x86. Es útil comenzar aquí para transponer matrices cuadradas más grandes, como 8x8 o 16x16.

_MM_TRANSPOSE4_PS(r0, r1, r2, r3)se implementa de manera diferente por diferentes compiladores. GCC e ICC (no he verificado Clang) usan unpcklps, unpckhps, unpcklpd, unpckhpdmientras que MSVC solo usa shufps. De hecho, podemos combinar estos dos enfoques juntos de esta manera.

t0 = _mm_unpacklo_ps(r0, r1);
t1 = _mm_unpackhi_ps(r0, r1);
t2 = _mm_unpacklo_ps(r2, r3);
t3 = _mm_unpackhi_ps(r2, r3);

r0 = _mm_shuffle_ps(t0,t2, 0x44);
r1 = _mm_shuffle_ps(t0,t2, 0xEE);
r2 = _mm_shuffle_ps(t1,t3, 0x44);
r3 = _mm_shuffle_ps(t1,t3, 0xEE);

Una observación interesante es que dos aleatorios se pueden convertir en uno aleatorio y dos combinaciones (SSE4.1) de esta manera.

t0 = _mm_unpacklo_ps(r0, r1);
t1 = _mm_unpackhi_ps(r0, r1);
t2 = _mm_unpacklo_ps(r2, r3);
t3 = _mm_unpackhi_ps(r2, r3);

v  = _mm_shuffle_ps(t0,t2, 0x4E);
r0 = _mm_blend_ps(t0,v, 0xC);
r1 = _mm_blend_ps(t2,v, 0x3);
v  = _mm_shuffle_ps(t1,t3, 0x4E);
r2 = _mm_blend_ps(t1,v, 0xC);
r3 = _mm_blend_ps(t3,v, 0x3);

Esto convirtió efectivamente 4 barajas en 2 barajas y 4 combinaciones. Esto usa 2 instrucciones más que la implementación de GCC, ICC y MSVC. La ventaja es que reduce la presión del puerto, lo que puede tener un beneficio en algunas circunstancias. Actualmente, todos los barajados y desempaquetados pueden ir solo a un puerto en particular, mientras que las mezclas pueden ir a cualquiera de dos puertos diferentes.

Intenté usar 8 aleatorios como MSVC y convertirlo en 4 aleatorios + 8 combinaciones, pero no funcionó. Todavía tenía que usar 4 desempaquetadores.

Usé esta misma técnica para una transposición flotante de 8x8 (ver hacia el final de esa respuesta). https://stackoverflow.com/a/25627536/2542702 . En esa respuesta, todavía tuve que usar 8 desempaquetados, pero logré convertir los 8 aleatorios en 4 aleatorios y 8 combinaciones.

Para los enteros de 32 bits, no hay nada parecido shufps(excepto para las mezclas de 128 bits con AVX512), por lo que solo se puede implementar con descomprimidos que no creo que se puedan convertir en mezclas (de manera eficiente). Con AVX512 vshufi32x4actúa de manera efectiva, shufpsexcepto para carriles de 128 bits de 4 enteros en lugar de flotantes de 32 bits, por lo que esta misma técnica podría usarse vshufi32x4en algunos casos. Con Knights Landing, las mezclas son cuatro veces más lentas (rendimiento) que las mezclas.

Bosón Z
fuente
1
Puede usar shufpsen datos enteros. Si está barajando mucho, podría valer la pena hacerlo todo en el dominio FP para shufps+ blendps, especialmente si no tiene vpblendddisponible el AVX2 igualmente eficiente . Además, en el hardware de la familia Intel SnB, no hay un retraso de derivación adicional para usar shufpsentre instrucciones enteras como paddd. (Hay un retraso de derivación para mezclar blendpscon paddd, de acuerdo con las pruebas de SNB Agner Fog, sin embargo.)
Peter Cordes
@PeterCordes, necesito revisar los cambios de dominio nuevamente. ¿Hay alguna tabla (tal vez una respuesta en SO) que resuma la penalización por cambio de dominio para Core2-Skylake? En cualquier caso, he pensado más en esto. Ahora veo por qué wim y tú seguían mencionando vinsertf64x4en mi respuesta de transposición de 16x16 en lugar de vinserti64x4. Si estoy leyendo y luego escribiendo la matriz, ciertamente no importa si uso el dominio de punto flotante o el dominio entero, ya que la transposición es solo datos en movimiento.
Bosón Z
1
Las tablas de Agner enumeran dominios por instrucción para Core2 y Nehalem (y creo que AMD), pero no la familia SnB. La guía de microarquía de Agner solo tiene un párrafo que dice que se reduce a 1c y, a menudo, 0 en SnB, con algunos ejemplos. Creo que el manual de optimización de Intel tiene una tabla, pero no he intentado asimilarla, así que no recuerdo cuántos detalles tiene. Recuerdo que no era totalmente obvio en qué categoría estaría una instrucción determinada.
Peter Cordes
Incluso si no solo está escribiendo en la memoria, es solo 1 reloj adicional para toda la transposición. El retardo adicional para cada operando puede ocurrir en paralelo (o de forma escalonada) cuando el consumidor de la transposición comienza a leer registros escritos por mezclas o mezclas. La ejecución fuera de orden permite que comiencen las primeras FMA o lo que sea mientras terminan las últimas barajas, pero no hay una cadena de demoras de dypass, solo una extra como máximo.
Peter Cordes
1
Nicw responde! El manual de optimización de arquitecturas de intel 64-ia-32, tabla 2-3, enumera los retrasos de omisión para Skylake, tal vez eso sea de su interés. La tabla 2-8 para Haswell se ve bastante diferente.
wim
1

Considere cada fila como una columna y cada columna como una fila ... use j, i en lugar de i, j

demostración: http://ideone.com/lvsxKZ

#include <iostream> 
using namespace std;

int main ()
{
    char A [3][3] =
    {
        { 'a', 'b', 'c' },
        { 'd', 'e', 'f' },
        { 'g', 'h', 'i' }
    };

    cout << "A = " << endl << endl;

    // print matrix A
    for (int i=0; i<3; i++)
    {
        for (int j=0; j<3; j++) cout << A[i][j];
        cout << endl;
    }

    cout << endl << "A transpose = " << endl << endl;

    // print A transpose
    for (int i=0; i<3; i++)
    {
        for (int j=0; j<3; j++) cout << A[j][i];
        cout << endl;
    }

    return 0;
}
Khaled.K
fuente
1

transposición sin gastos generales (clase no completa):

class Matrix{
   double *data; //suppose this will point to data
   double _get1(int i, int j){return data[i*M+j];} //used to access normally
   double _get2(int i, int j){return data[j*N+i];} //used when transposed

   public:
   int M, N; //dimensions
   double (*get_p)(int, int); //functor to access elements  
   Matrix(int _M,int _N):M(_M), N(_N){
     //allocate data
     get_p=&Matrix::_get1; // initialised with normal access 
     }

   double get(int i, int j){
     //there should be a way to directly use get_p to call. but i think even this
     //doesnt incur overhead because it is inline and the compiler should be intelligent
     //enough to remove the extra call
     return (this->*get_p)(i,j);
    }
   void transpose(){ //twice transpose gives the original
     if(get_p==&Matrix::get1) get_p=&Matrix::_get2;
     else get_p==&Matrix::_get1; 
     swap(M,N);
     }
}

se puede utilizar así:

Matrix M(100,200);
double x=M.get(17,45);
M.transpose();
x=M.get(17,45); // = original M(45,17)

por supuesto, no me molesté con la gestión de la memoria aquí, que es un tema crucial pero diferente.

Reza Baram
fuente
4
Tiene una sobrecarga de su puntero de función que debe seguirse para cada acceso de elemento.
user877329
1

Si conocemos el tamaño de las matrices antes, podríamos utilizar la unión como ayuda. Me gusta esto-

#include <bits/stdc++.h>
using namespace std;

union ua{
    int arr[2][3];
    int brr[3][2];
};

int main() {
    union ua uav;
    int karr[2][3] = {{1,2,3},{4,5,6}};
    memcpy(uav.arr,karr,sizeof(karr));
    for (int i=0;i<3;i++)
    {
        for (int j=0;j<2;j++)
            cout<<uav.brr[i][j]<<" ";
        cout<<'\n';
    }

    return 0;
}
Sandeep KV
fuente
Soy nuevo en C / C ++, pero esto parece genial. Debido a que union usa la ubicación de memoria compartida para sus miembros, puede leer esa memoria de manera diferente. Por lo tanto, obtiene una matriz transpuesta sin hacer una nueva asignación de matriz. Estoy en lo cierto?
Doğuş
1
template <class T>
void transpose( const 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];
        }
    }
} 
Raquel gallen
fuente
1
Prefiero pensar que sería más rápido si intercambia los dos bucles, debido a una menor penalización por error de caché al escribir que al leer.
phoeagon
5
Esto solo funciona para una matriz cuadrada. ¡Una matriz rectangular es un problema completamente diferente!
NealB
2
La pregunta pide la forma más rápida. Esta es solo una forma. ¿Qué te hace pensar que es rápido, y mucho menos más rápido? Para matrices grandes, esto destruirá el caché y tendrá un rendimiento terrible.
Eric Postpischil
1
@NealB: ¿Cómo se imagina eso?
Eric Postpischil
@EricPostpischil El OP pregunta por una matriz relativamente grande, así que supongo que querían hacerlo "en su lugar" para evitar asignar el doble de memoria. Cuando se hace esto, la dirección base de las matrices de origen y destino es la misma. La transposición cambiando los índices de fila y columna solo funcionará para matrices cuadradas. Hay métodos para hacerlo bien en matrices rectangulares, pero son algo más complejos.
NealB
0

Las bibliotecas de álgebra lineal modernas incluyen versiones optimizadas de las operaciones más comunes. Muchos de ellos incluyen distribución dinámica de CPU, que elige la mejor implementación para el hardware en el momento de ejecución del programa (sin comprometer la portabilidad).

Esta es comúnmente una mejor alternativa para realizar la optimización manual de sus functinos a través de funciones intrínsecas de extensiones vectoriales. Este último vinculará su implementación a un proveedor y modelo de hardware en particular: si decide cambiar a un proveedor diferente (por ejemplo, Power, ARM) oa una extensión de vector más reciente (por ejemplo, AVX512), deberá volver a implementarlo nuevamente para aprovecharlos al máximo.

La transposición MKL, por ejemplo, incluye la función de extensiones BLAS imatcopy. También puede encontrarlo en otras implementaciones como OpenBLAS:

#include <mkl.h>

void transpose( float* a, int n, int m ) {
    const char row_major = 'R';
    const char transpose = 'T';
    const float alpha = 1.0f;
    mkl_simatcopy (row_major, transpose, n, m, alpha, a, n, n);
}

Para un proyecto de C ++, puede hacer uso de Armadillo C ++:

#include <armadillo>

void transpose( arma::mat &matrix ) {
    arma::inplace_trans(matrix);
}
Jorge Bellon
fuente
0

intel mkl sugiere matrices de transposición / copia en el lugar y fuera de lugar. aquí está el enlace a la documentación . Recomendaría probar la implementación fuera de lugar ya que diez más rápido en el lugar y en la documentación de la última versión de mkl contiene algunos errores.

Gennady.F
fuente
-1

Creo que la forma más rápida no debería tomar más alto que O (n ^ 2) también de esta manera puede usar solo el espacio O (1):
la forma de hacerlo es intercambiar en pares porque cuando transpone una matriz, entonces lo que hacer es: M [i] [j] = M [j] [i], así que almacene M [i] [j] en temp, luego M [i] [j] = M [j] [i], y el último paso: M [j] [i] = temp. esto podría hacerse con una sola pasada, por lo que debería tomar O (n ^ 2)

Fayez Abdlrazaq Deab
fuente
2
M [i] [j] = M [j] [i] solo funcionará si fuera una matriz cuadrada; de lo contrario, arrojaría una excepción de índice.
Antony Thomas
-6

mi respuesta está transpuesta de matriz 3x3

 #include<iostream.h>

#include<math.h>


main()
{
int a[3][3];
int b[3];
cout<<"You must give us an array 3x3 and then we will give you Transposed it "<<endl;
for(int i=0;i<3;i++)
{
    for(int j=0;j<3;j++)
{
cout<<"Enter a["<<i<<"]["<<j<<"]: ";

cin>>a[i][j];

}

}
cout<<"Matrix you entered is :"<<endl;

 for (int e = 0 ; e < 3 ; e++ )

{
    for ( int f = 0 ; f < 3 ; f++ )

        cout << a[e][f] << "\t";


    cout << endl;

    }

 cout<<"\nTransposed of matrix you entered is :"<<endl;
 for (int c = 0 ; c < 3 ; c++ )
{
    for ( int d = 0 ; d < 3 ; d++ )
        cout << a[d][c] << "\t";

    cout << endl;
    }

return 0;
}
ángel
fuente