Saya punya matrix (relatif besar) yang saya butuhkan untuk merefleksikan. Misalnya menganggap bahwa saya adalah matriks
a b c d e f
g h i j k l
m n o p q r
Saya ingin hasilnya menjadi sebagai berikut:
a g m
b h n
c I o
d j p
e k q
f l r
Apa cara tercepat untuk melakukan ini?
Ini adalah pertanyaan yang bagus. Ada banyak alasan anda ingin benar-benar transpos matriks dalam memori bukan hanya swap koordinat, misalnya dalam matriks perkalian dan Gaussian mengolesi.
Pertama ijinkan saya menyebutkan salah satu fungsi yang saya gunakan untuk transpose (EDIT: silahkan lihat akhir dari jawaban saya di mana saya menemukan sebuah solusi jauh lebih cepat)
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];
}
}
Sekarang mari's melihat mengapa transpose berguna. Mempertimbangkan perkalian matriks C = A*B. Kita bisa melakukannya dengan cara ini.
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;
}
}
Dengan cara itu, bagaimanapun, akan memiliki banyak cache misses. Jauh lebih cepat solusinya adalah untuk mengambil transpos dari B pertama
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);
Perkalian matriks adalah O(n^3) dan transpose adalah O(n^2), sehingga mengambil transpose harus memiliki efek yang dapat diabaikan pada perhitungan waktu (untuk besar n
). Dalam matriks perkalian loop ubin bahkan lebih efektif daripada mengambil transpose tapi yang's jauh lebih rumit.
Aku berharap aku tahu cara yang lebih cepat untuk melakukan transpose (Edit: saya menemukan solusi yang lebih cepat, melihat akhir dari jawaban saya). Ketika Haswell/AVX2 keluar dalam beberapa minggu ini akan memiliki berkumpul fungsi. Saya don't tahu apakah itu akan membantu dalam kasus ini tapi aku bisa gambar pertemuan kolom dan menulis satu baris. Mungkin itu akan membuat transpose yang tidak perlu.
Untuk Gaussian mengolesi apa yang anda lakukan adalah mengolesi horizontal dan kemudian smear secara vertikal. Tapi mengolesi secara vertikal memiliki cache masalah jadi apa yang anda lakukan adalah
Smear image horizontally
transpose output
Smear output horizontally
transpose output
Berikut ini adalah kertas oleh Intel menjelaskan bahwa http://software.intel.com/en-us/articles/iir-gaussian-blur-filter-implementation-using-intel-advanced-vector-extensions
Terakhir, apa yang saya benar-benar melakukan dalam perkalian matriks (dan Gaussian mengolesi) tidak mengambil persis transpose, tapi mengambil transpose dalam lebar tertentu ukuran vektor (misalnya 4 atau 8 untuk SSE/AVX). Berikut adalah fungsi yang saya gunakan
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:
Aku mencoba beberapa fungsi untuk menemukan cara tercepat transpose untuk matriks besar. Pada akhirnya hasil tercepat adalah dengan menggunakan loop memblokir dengan block_size=16
(Edit: saya menemukan solusi yang lebih cepat menggunakan SSE dan loop blocking - lihat di bawah). Kode ini bekerja untuk setiap matriks NxM (yaitu matriks tidak harus persegi).
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);
}
}
}
Nilai lda
dan ldb
adalah lebar dari matriks. Ini harus kelipatan dari ukuran blok. Untuk menemukan nilai-nilai dan mengalokasikan memori untuk misalnya 3000x1001 matrix saya melakukan sesuatu seperti ini
#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);
Untuk 3000x1001 ini kembali ldb = 3008
dan lda = 1008
Edit:
Saya menemukan bahkan solusi cepat menggunakan SSE intrinsik:
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);
}
}
}
}
}
Beberapa rincian tentang transposing 4x4 meter mengapung (saya akan membahas 32-bit integer kemudian) matriks dengan perangkat keras x86. It's membantu untuk mulai di sini dalam rangka untuk transpose matriks persegi yang lebih besar seperti 8x8 atau 16 x 16.
_MM_TRANSPOSE4_PS(r0, r1, r2, r3)
yang diterapkan berbeda-beda oleh berbagai compiler. GCC dan ICC (aku belum diperiksa Dentang) menggunakan unpcklps, unpckhps, unpcklpd, unpckhpd
sedangkan MSVC hanya menggunakan shufps
. Kita benar-benar dapat menggabungkan dua pendekatan bersama-sama seperti ini.
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);
Salah satu pengamatan yang menarik adalah bahwa dua mengocok dapat dikonversi ke satu shuffle dan dua campuran (SSE4.1) seperti ini.
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);
Ini secara efektif diubah 4 mengocok menjadi 2 mengaduk-aduk dan 4 campuran. Ini menggunakan 2 petunjuk lebih lanjut dari pelaksanaan GCC, ICC, dan MSVC. Keuntungan adalah bahwa itu mengurangi tekanan pelabuhan yang mungkin memiliki manfaat dalam beberapa keadaan. Saat ini semua mengocok dan membongkar bisa pergi hanya untuk satu port tertentu sedangkan campuran dapat pergi ke salah satu dari dua port yang berbeda.
Saya mencoba menggunakan 8 mengaduk-aduk seperti MSVC dan mengkonversi ke 4 mengocok + 8 menyatu tapi itu tidak bekerja. Saya masih harus menggunakan 4 membongkar.
Saya menggunakan teknik yang sama untuk 8x8 mengapung transpose (lihat menjelang akhir dari yang menjawab). https://stackoverflow.com/a/25627536/2542702. Dalam jawaban itu saya masih harus menggunakan 8 membongkar tapi aku manged untuk mengkonversi 8 mengocok menjadi 4 mengocok dan 8 menyatu.
Untuk 32-bit bilangan bulat tidak ada yang seperti shufps
(kecuali untuk 128-bit mengocok dengan AVX512) sehingga hanya dapat dilaksanakan dengan membongkar yang saya don't pikir dapat mengkonversi ke menyatu (efisien). Dengan AVX512 vshufi32x4
bertindak secara efektif seperti shufps
kecuali untuk 128-bit jalur 4 bukan bilangan bulat 32-bit mengapung sehingga teknik yang sama ini mungkin mungkin dengan vshufi32x4
dalam beberapa kasus. Dengan Knights Landing mengocok empat kali lebih lambat (throughput) dari campuran.
transposing tanpa overhead (kelas tidak lengkap):
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);
}
}
dapat digunakan seperti ini:
Matrix M(100,200);
double x=M.get(17,45);
M.transpose();
x=M.get(17,45); // = original M(45,17)
tentu saja aku didn't repot-repot dengan manajemen memori di sini, yang penting tapi topik yang berbeda.
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];
}
}
}
Pertimbangkan setiap baris sebagai kolom, dan setiap kolom sebagai turut .. menggunakan j,aku bukan aku,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;
}
intel mkl menunjukkan di tempat dan out-of-tempat transposisi/menyalin matriks. berikut adalah link ke dokumentasi. Saya akan merekomendasikan mencoba keluar dari tempat pelaksanaan serta lebih cepat sepuluh di tempat dan menjadi dokumentasi dari versi terbaru dari mkl mengandung beberapa kesalahan.
Jika ukuran dari array yang diketahui sebelumnya maka kita bisa menggunakan union untuk membantu kami. Seperti ini-
``
menggunakan namespace std;
uni ua{ int arr[2][3]; int brr[3][2]; };
int main() { uni 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; } ``
Modern aljabar linear perpustakaan mencakup dioptimalkan versi yang paling umum operasi. Banyak dari mereka termasuk dynamic CPU pengiriman, yang memilih implementasi terbaik untuk perangkat keras pada program waktu pelaksanaan (tanpa mengorbankan portabilitas).
Hal ini sering alternatif yang lebih baik untuk melakukan optimasi manual anda functinos melalui vector extensions fungsi intrinsik. Yang terakhir akan mengikat implementasi tertentu vendor hardware dan model: jika anda memutuskan untuk swap untuk vendor yang berbeda (misalnya Kekuatan, LENGAN) atau yang lebih baru vector extensions (misalnya AVX512), anda akan perlu untuk menerapkan kembali lagi untuk mendapatkan yang paling dari mereka.
MKL transposisi, misalnya, termasuk BLAS ekstensi fungsi imatcopy
. Anda dapat menemukannya dalam implementasi lainnya seperti OpenBLAS juga:
``
batal transpose( float* a, int n, int m ) { const char row_major = 'R'; const char transpose = 'T'; const float alpha = 1.0 f; mkl_simatcopy (row_major, transpose, n, m, alpha, a, n, n); } ``
Untuk C++ proyek, anda dapat menggunakan Armadillo C++: ``
batal transpose( arma::mat &matriks ) { arma::inplace_trans(matrix); } ``
Saya pikir yang paling cepat cara tidak harus mengambil lebih tinggi dari O(n^2) juga dengan cara ini anda dapat menggunakan hanya O(1) ruang : cara untuk melakukannya adalah untuk swap di pasang karena ketika anda transpose matriks maka apa yang anda lakukan adalah: M[i][j]=M[j][i] , sehingga toko M[i][j] di temp, kemudian M[i][j]=M[j][i],dan langkah terakhir : M[j][i]=temp. hal ini bisa dilakukan oleh satu lulus sehingga harus mengambil O(n^2)
jawaban saya adalah dialihkan dari 3x3 matrix
#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;
}