Какой самый быстрый способ транспонировать матрицу в C++?
У меня есть матрица (относительно большая), которую мне нужно транспонировать. Например предположим, что моя матрица
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
Какой самый быстрый способ сделать это?
12 ответов
Это хороший вопрос. Есть много причин, по которым вы захотите переставить матрицу в памяти, а не просто поменять координаты, например, при умножении матрицы и смазывании по Гауссу.
Сначала позвольте мне перечислить одну из функций, которые я использую для транспонирования (РЕДАКТИРОВАТЬ: см. В конце моего ответа, где я нашел гораздо более быстрое решение)
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). В матричном умножении циклическое разбиение еще более эффективно, чем использование транспонирования, но это намного сложнее.
Хотелось бы, чтобы я знал более быстрый способ сделать транспонирование (Edit: я нашел более быстрое решение, см. Конец моего ответа). Когда через несколько недель выйдет 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 ширина матрицы. Они должны быть кратны размеру блока. Чтобы найти значения и выделить память, например, для матрицы 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
Редактировать:
Я нашел еще более быстрое решение с использованием встроенных функций 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 (см. В конце этого ответа). /questions/388340/transponirovat-poplavok-8x8-ispolzuya-avxavx2/388363#388363. В этом ответе мне все еще пришлось использовать 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;
}
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];
}
}
}
Рассматривайте каждую строку как столбец, а каждый столбец - как строку. Используйте 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)
конечно, я не беспокоился об управлении памятью здесь, что является важной, но другой темой.
Intel mkl предлагает матрицы транспонирования / копирования на месте и вне места. вот ссылка на документацию. Я бы рекомендовал попробовать неуместную реализацию, так как более быстрая десятка на месте и в документации последней версии mkl есть некоторые ошибки.
Современные библиотеки линейной алгебры включают оптимизированные версии наиболее распространенных операций. Многие из них включают динамическую диспетчеризацию ЦП, которая выбирает лучшую реализацию для оборудования во время выполнения программы (без ущерба для переносимости).
Обычно это лучшая альтернатива ручной оптимизации ваших функций с помощью встроенных функций векторных расширений. Последнее свяжет вашу реализацию с конкретным поставщиком оборудования и моделью: если вы решите перейти на другого поставщика (например, 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);
}
Самая быстрая транспозиция — это та, которая останется в кеше для следующей операции (которая будет ее использовать).
Например, не транспонируйте все сразу. Транспонируйте только одну подматрицу. Затем используйте его в части следующего алгоритма, которому нужны транспонированные данные. Затем транспонируйте следующую подматрицу. Потом посчитай. Затем транспонируйте другую подматрицу. Повторяйте, пока вся матрица не будет транспонирована. Таким образом, данные сохраняются в кеше.
Если вы сразу полностью транспонируете 128-мегабайтную матрицу на процессор с 2-мегабайтным кешем, то в конце операции в кеше останутся только последние биты матрицы. Тогда вам лучше начать умножать матрицы из последних битов, чтобы использовать эти горячие данные размером 2 МБ.
Но когда вы разделите работу на более мелкие части, например, умножив подматрицы, вы можете просто выполнить ленивую транспозицию следующим образом:
multiply:
for all sub_matrices in mat1 row
for all sub_matrices in mat2 column
select sub_matrix1
select sub_matrix2
if sub_mat2 is not transposed
transpose sub_mat2
multiply sub_mat1 and sub_mat2 <---- data in cache!
accumulate result
Преимущества:
- Пропускная способность кэша L1/L2 используется для следующей операции.
- задержка транспозиции скрыта за вычислением следующей операции
- работает с небольшим кешем, всего 64 КБ, зависит от размера чанка
Я думаю, что самый быстрый способ не должен брать больше, чем O(n^2), и таким образом вы можете использовать только O(1) пробел:
способ сделать это - поменяться парами, потому что когда вы перемещаете матрицу, то вы делаете следующее: M[i] [j]= M [j] [i], поэтому сохраняйте M [i] [j] в temp, тогда M [i] [j]= M [j] [i], и последний шаг: M[j] [i]= темп. это может быть сделано за один проход, поэтому это должно занять O(n^2)
Мой ответ транспонирован из матрицы 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;
}