Каков самый быстрый способ транспонирования матрицы в C++?

Question

У меня есть матрица (относительно большая), которую мне нужно транспонировать. Например предположим, что моя матрица

a b c d e f 
g h i j k l 
m n o p q r 

Я хочу, чтобы результат будет следующим:

a g m 
b h n 
c I o 
d j p 
e k q 
f l r 

Какой самый быстрый способ сделать это?

Answer 1

Это хороший вопрос. Существует много причин, по которым вы хотели бы фактически перенести матрицу в память, а не просто менять координаты, например. в матричном умножении и гауссовском размытии.

Прежде всего позвольте мне перечислить одну из функций, которые я использую для транспонирования (EDIT: см конец моего ответа, где я нашел гораздо более быстрое разрешение)

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

Теперь давайте посмотрим, почему транспонированная Полезно. Рассмотрим матричное умножение C = A * 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*l+j]; 
     } 
     C[K*i + j] = tmp; 
    } 
} 

Этот способ, однако, будет иметь много недостатков в кэше. Гораздо быстрее, решение взять транспонирование В первой

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); 

умножения матриц О (п^3) и транспонирования O (N^2), поэтому принимать транспонирование должны иметь незначительное влияние на время вычисления (для больших n). В матричном умножении петля еще более эффективна, чем перенос, но это намного сложнее.

Мне жаль, что я не знал более быстрый способ сделать транспонирование (Редактировать: я нашел более быстрое решение, см. Конец моего ответа). Когда Haswell/AVX2 выйдет через несколько недель, у него будет функция сбора. Я не знаю, будет ли это полезно в этом случае, но я мог бы собрать столбец и написать строку. Возможно, это сделает ненужным транспозицию.

Для гауссова размытия, что вы делаете, размазывается горизонтально, а затем размазывается вертикально. Но размазывание по вертикали имеет проблемы кэша, так что вы делаете

Smear image horizontally 
transpose output 
Smear output horizontally 
transpose output 

Вот документ от Intel, объясняя, что http://software.intel.com/en-us/articles/iir-gaussian-blur-filter-implementation-using-intel-advanced-vector-extensions

Наконец, то, что я на самом деле делать в матричное умножение (и в гауссовой размытию) не возьмите именно транспонирование, но возьмите транспонирование по ширине определенного векторного размера (например, 4 или 8 для SSE/AVX). Вот функция, я использую

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:

Я попробовал несколько функций, чтобы найти самый быстрый транспонировать для больших матриц.В итоге самым быстрым результатом является использование блокировки цикла с помощью block_size=16 (Редактирование: я нашел более быстрое решение с использованием SSE и блокировки циклов - см. Ниже). Этот код работает для любой матрицы NxM (т. Е. Матрица не должна быть квадратной).

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

Значения lda и ldb являются ширина матрицы. Они должны быть кратными размеру блока. Чтобы найти значения и выделить память, например. 3000x1001 матрица я что-то вроде этого

#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); 

Для 3000x1001 это возвращает ldb = 3008 и lda = 1008

Edit:

я нашел еще более быстрое решение с использованием 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); 
       } 
      } 
     } 
    } 
} 

Answer 2

Это будет зависеть от вашего приложения, но в целом самый быстрый способ транспонировать матрицу будет состоять в том, чтобы инвертировать ваши координаты, когда вы смотрите вверх, тогда вам не нужно фактически перемещать какие-либо данные.

Answer 3

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

Answer 4

Рассмотрим каждую строку как столбец и каждый столбец в виде строки .. используйте j, i inst Свинец из I, J

демо: 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; 
} 

Answer 5

Я думаю, что самый быстрый способ не принимать больше, чем O (N^2) и, таким образом, вы можете использовать только O (1) пространство:
способ сделать это, чтобы поменять местами, потому что, когда вы транспонируете матрицу, тогда вы делаете следующее: M [i] [j] = M [j] [i], поэтому сохраните M [i] [j] в temp, то M [i] [j] = M [j] [i], и последний шаг: M [j] [i] = temp. это может быть сделано за один проход, так что следует принимать O (N^2)

Answer 6

мой ответ транспонируется в виде матрицы 3х3

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

Answer 7

Транспонирование без каких-либо накладных расходов (класс не полный):

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

можно использовать так:

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

я, конечно, не заморачиваться с управлением памятью здесь, что очень важно, но другая тема.

Answer 8

Некоторые сведения о транспонировании квадратного поплавка 4x4 (я расскажу позже о 32-битном целочисленном) с использованием оборудования x86. Здесь полезно начать транспонировать большие квадратные матрицы, такие как 8x8 или 16x16.

_MM_TRANSPOSE4_PS(r0, r1, r2, r3) осуществляется по-разному различными компиляторами. GCC и ICC (я не проверял Clang) используют unpcklps, unpckhps, unpcklpd, unpckhpd, тогда как MSVC использует только shufps. Мы можем объединить эти два подхода вместе.

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); 

Одно интересное наблюдение состоит в том, что два перетасовки могут быть преобразованы в один перетасовки и двух смесей (SSE4.1), как это.

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); 

Это эффективно преобразованное 4 тасования в 2 тасования и 4 смеси. Это использует еще 2 инструкции, чем реализация GCC, ICC и MSVC. Преимущество заключается в том, что он снижает давление порта, которое может принести пользу в некоторых случаях. В настоящее время все перетасовки и распаковки могут идти только на один конкретный порт, тогда как смеси могут перейти в любой из двух разных портов.

Я попытался использовать 8 тасов, таких как MSVC, и преобразовал их в 4 тасования + 8 смесей, но это не сработало. Мне все еще пришлось использовать 4 распаковки.

Я использовал эту же методику для транспортера с поплавком 8x8 (см. Конец этого ответа). https://stackoverflow.com/a/25627536/2542702. В этом ответе мне все еще приходилось использовать 8 распаковщиков, но мне удалось преобразовать 8 тасов в 4 тасования и 8 смесей.

Для 32-разрядных целых чисел нет ничего, кроме shufps (за исключением 128-битных перетасовки с AVX512), поэтому его можно реализовать только с помощью распаковки, которые, как я думаю, не могут быть конвертированы в смеси (эффективно). С AVX512 vshufi32x4 действует эффективно, как shufps, за исключением 128-битных дорожек из 4 целых чисел вместо 32-битных поплавков, поэтому в некоторых случаях этот метод может быть, возможно, с vshufi32x4. С рыцарями Посадка в случайном порядке в четыре раза медленнее (пропускная способность), чем смеси.