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

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

Сначала позвольте мне перечислить одну из функций, которые я использую для транспонирования (РЕДАКТИРОВАТЬ: пожалуйста, посмотрите конец моего ответа, где я нашел гораздо более быстрое решение)

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

Таким образом, однако, будет много промахов в кеше. Гораздо более быстрое решение - сначала перенести букву 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);

Умножение матрицы - O (n ^ 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];
    }
}

РЕДАКТИРОВАТЬ:

Я попробовал несколько функций, чтобы найти наиболее быстрое транспонирование для больших матриц. В конце концов, самым быстрым результатом является использование блокировки цикла с 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 - это ширина матрицы. Они должны быть кратны размеру блока. Чтобы найти значения и выделить память, например, матрица 3000х1001 делаю примерно так

#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

Изменить:

Я нашел еще более быстрое решение с использованием встроенных функций 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);
                }
            }
        }
    }
}

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

Некоторые подробности о транспонировании квадратных матриц 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. С Knights Landing перемешивание в четыре раза медленнее (пропускная способность), чем смешивание.

Если размер массивов известен заранее, мы могли бы использовать объединение для нашей помощи. Нравится-

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

Рассматривайте каждую строку как столбец, а каждый столбец как строку ... используйте j, i вместо 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;
}

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

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)

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

Современные библиотеки линейной алгебры включают оптимизированные версии наиболее распространенных операций. Многие из них включают динамическую диспетчеризацию ЦП, которая выбирает лучшую реализацию для оборудования во время выполнения программы (без ущерба для переносимости).

Обычно это лучшая альтернатива ручной оптимизации ваших функций с помощью встроенных функций векторных расширений. Последнее свяжет вашу реализацию с конкретным поставщиком оборудования и моделью: если вы решите переключиться на другого поставщика (например, Power, ARM) или на более новые векторные расширения (например, AVX512), вам нужно будет повторно реализовать его снова, чтобы получить от них максимум.

Например, транспонирование MKL включает функцию расширений BLAS imatcopy. Вы также можете найти его в других реализациях, таких как 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);
}

Для проекта C ++ вы можете использовать Armadillo C ++:

#include <armadillo>

void transpose( arma::mat &matrix ) {
    arma::inplace_trans(matrix);
}

intel mkl предлагает матрицы транспонирования / копирования на месте и не на месте. вот ссылка на документацию. Я бы рекомендовал попробовать неуместную реализацию, так как более быстрая десятка на месте и в документации последней версии mkl есть некоторые ошибки.

Я думаю, что самый быстрый способ не должен превышать 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)

мой ответ транспонирован из матрицы 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;
}