У меня есть матрица (относительно большая), которую мне нужно транспонировать. Например, предположим, что моя матрица имеет вид
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
Какой самый быстрый способ сделать это?
Это хороший вопрос. Есть много причин, по которым вы захотите действительно транспонировать матрицу в памяти, а не просто поменять координаты местами, например, при матричном умножении и гауссовом размазывании.
Сначала позвольте мне перечислить одну из функций, которые я использую для транспонирования (ПРИМЕЧАНИЕ: пожалуйста, смотрите конец моего ответа, где я нашел гораздо более быстрое решение)
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, в нем будет функция gather. Я не знаю, будет ли она полезна в данном случае, но я могу представить, как я собираю столбец и выписываю строку. Возможно, это сделает транспонирование ненужным.
Для гауссова размазывания вы делаете горизонтальное размазывание, а затем вертикальное размазывание. Но при вертикальном размазывании возникает проблема с кэшем, поэтому нужно сделать следующее.
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
Редактирование:
Я нашел еще более быстрое решение, используя SSE intrinsics:
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);
}
}
}
}
}
Это зависит от вашего приложения, но в целом самым быстрым способом транспонирования матрицы будет инвертирование координат при поиске, тогда вам не придется перемещать данные.
Некоторые подробности о переносе 4х4 квадратных поплавок (я буду обсуждать 32-разрядное целое число позже) матриц с x86-оборудовании. Это's полезн для начала здесь для того, чтобы перенести больше площади матриц 8х8 или 16х16.
_MM_TRANSPOSE4_PS(Р0, Р1, Р2, Р3)
реализуется по-разному в разных компиляторах. На GCC и ICC (я не проверял лязг) 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 инструкций, чем реализация ССЗ, МУС, и MSVC. Преимущество в том, что он уменьшает давление, которое может иметь преимущества в некоторых случаях. В настоящее время все перемешивает и распаковывает можете пойти только на один конкретный порт, а смеси может перейти к любому из двух разных портов.
Я попытался с помощью 8 тасует, как MSVC и преобразования в 4 тасует + 8 блендов, но он не работал. Я все еще должен был использовать 4 распаковывает.
Я использовал этот же метод для 8х8 плавающий транспонировать (см. В конце ответа). https://stackoverflow.com/a/25627536/2542702. В том, что ответа я все равно пришлось использовать 8 распаковывает но я manged, чтобы преобразовать 8 тасует в 4 тасует и 8 блендов.
Для 32-разрядных целых чисел нет ничего похожего на shufps (за исключением 128-битное тасует с AVX512), поэтому она может быть осуществлена только с распаковывает, которые я не'т думаю, что можно преобразовать в смеси (эффективно). Эффективно действует AVX512
vshufi32x4 " как " shufps, кроме 128-битное полосы из 4 чисел вместо 32-разрядных плавает, так этот же метод может быть возможно с vshufi32x4 в некоторых случаях. С рыцарями посадки тасует в четыре раза медленнее (пропускная способность), чем смеси.
транспонирование без каких-либо накладных расходов (класс не полный):
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)
конечно, я этого'т возиться с управление памяти вот, что имеет решающее значение, но разные темы.
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;
}
Интел мкл указывает на место и место переноса/копирования матриц. вот ссылка на документацию. Я бы рекомендовал пробовать из места внедрения как быстрее десяти на месте и в документации последняя версия мкл содержит некоторые ошибки.
Если размер массива известен до того, как мы могли бы использовать союз с нашей помощью. Как это-
``
с помощью пространства имен std;
Союз уа{ инт модуль arr[2][3]; инт брр[3][2]; };
тап_п() { Союз уа БПЛА; инт Карр[2][3] = {{1,2,3},{4,5,6}}; функции memcpy(БПЛА.Арр Карр как sizeof(Карр)); для (int я=0;я&Л;3;я++) { для (Int J=0 и;ж<2;к++) соиь<<БПЛА.брр[я][Дж]<<" по себе "; соиь<<'\П'; }
возврат 0; } ``
Современные библиотеки линейной алгебры включают в себя оптимизированные версии из самых распространенных операций. Многие из них включают динамическая диспетчеризация процессора, который выбирает лучшее внедрение оборудования во время выполнения программы (без ущерба для портативности).
Обычно это лучшая альтернатива ручной оптимизации вашего functinos через вектор расширения встроенных функций. Последний будет связывать реализацию конкретного производителя оборудования и модели: если вы решите поменять поставщика (например, руку) или до новой векторных расширений (например, AVX512), вам придется повторно реализовать его снова, чтобы получить большинство из них.
МКЛ транспозиции, например, включает в себя imatcopy функция расширения Блас. Вы можете найти его в других реализациях, таких как OpenBLAS, а также:
пустота транспонировать( поплавок* а, инт Н, инт м ) { константный тип char row_major = 'Р'; константный тип char транспонирует = 'Т'; константный поплавок Альфа = 1.0 Ф; mkl_simatcopy (row_major, транспонировать, Н, М, АЛЬФА, а, н, н); } ``
Для проекта c++, вы можете воспользоваться броненосца на C++: ``
пустота транспонировать( АРМА::мат &матрица ) { АРМА::inplace_trans(матрицы); } ``
Я думаю, что самый быстрый способ не следует принимать выше, чем О(П^2) аналогичным способом можно использовать только O(1) по пространству : лучший способ сделать это, чтобы поменять в паре, потому что, когда вы транспонировать матрицу, то что вы делаете это: м[я][Дж]=м[Дж][я] , так что магазин, М[я][Дж] в темп, то m[я][Дж]=м[Дж][я],и последний шаг : м[Дж][я]=темп. это может быть сделано за один проход, так что его следует принимать за О(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;
}