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

Ответ 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;
    }
}

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

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

Ответ 2

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

Ответ 3

Некоторые подробности о транспонировании квадратного с плавающей запятой 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 распаковки.

Я использовал эту же технику для транспондера float 8x8 (см. конец этого ответа). fooobar.com/questions/4578/.... В этом ответе мне все еще приходилось использовать 8 распаковщиков, но мне удалось преобразовать 8 тасов в 4 тасования и 8 смесей.

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

Ответ 4

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

Ответ 5

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

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

Ответ 6

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

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)

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

Ответ 7

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

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

Ответ 8

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

Обычно это лучшая альтернатива выполнению ручной оптимизации ваших функций через встроенные функции векторных расширений. Последний будет привязывать вашу реализацию к конкретному поставщику и модели оборудования: если вы решите поменяться с другим поставщиком (например, 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++ вы можете использовать броненосец C++:

#include <armadillo>

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

Ответ 9

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

Ответ 10

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

Ответ 11

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