私は転置する必要がある行列(比較的大きい)を持っています。例えば、私の行列が次のようなものだとします。
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
一番手っ取り早い方法は何でしょうか?
これは良い質問ですね。 例えば、行列の乗算やガウススミアなど、座標を入れ替えるだけでなく、実際にメモリ上で行列を転置したい理由はたくさんあります。
まず、私が転置に使っている関数の1つを挙げてみましょう(編集:もっと速い解決策を見つけたので、回答の最後をご覧ください)。
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
以下は、それを説明するインテルの論文です。 http://software.intel.com/en-us/articles/iir-gaussian-blur-filter-implementation-using-intel-advanced-vector-extensions
最後に、行列の乗算で実際にやっていること(ガウススミアでも)は、正確に転置を取るのではなく、あるベクトルサイズ(SSE/AVXでは4や8など)の幅で転置を取ります。 以下は私が使っている関数です。
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</code> が返ってきます;
Edit:
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);
}
}
}
}
}
x86ハードウェアを使用した4x4スクエアフロート(後で32ビット整数について説明します)マトリックスの転置に関する詳細。 8x8や16x16などの大きな正方形の行列を転置するには、ここから始めると便利です。
_MM_TRANSPOSE4_PS(r0、r1、r2、r3)
は、コンパイラによって異なる方法で実装されます。 GCCとICC(私はClangをチェックしていません)は「unpcklps、unpckhps、unpcklpd、unpckhpd」を使用しますが、MSVCは「shufps」のみを使用します。 このように、これら2つのアプローチを実際に組み合わせることができます。
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);
興味深い観察の1つは、2つのシャッフルを1つのシャッフルと2つのブレンド(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つのブレンドに効果的に変換されました。 これには、GCC、ICC、およびMSVCの実装よりも2つ多くの命令を使用します。利点は、状況によってはメリットがあるポート圧力を軽減することです。 現在、すべてのシャッフルとアンパックは特定のポートにのみ移動できますが、ブレンドは2つの異なるポートのいずれかに移動できます。
MSVCのような8つのシャッフルを使用して、4つのシャッフル+ 8つのブレンドに変換してみましたが、機能しませんでした。 私はまだ4つのアンパックを使用する必要がありました。
これと同じ手法を8x8フロートトランスポーズに使用しました(その回答の最後を参照)。 https://stackoverflow.com/a/25627536/2542702。 その答えでは、まだ8つのアンパックを使用する必要がありましたが、8つのシャッフルを4つのシャッフルと8つのブレンドに変換するように管理しました。
32ビット整数の場合、「shufps」のようなものはありません(AVX512を使用した128ビットシャッフルを除く)。そのため、ブレンドに変換できない(効率的に)アンパックでのみ実装できます。 AVX512では、 vshufi32x4
は32ビットフロートの代わりに4整数の128ビットレーンを除いて shufps
のように効果的に作用するため、この同じ手法が vshufi32x4
の場合もあります。 ナイツランディングでは、シャッフルはブレンドの4倍遅い(スループット)。
オーバーヘッドなしで転置(クラスは完了していません):
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];
}
}
}
各行を列として、各列を行として考えます。 .. i、jの代わりにj、iを使用します。
#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は、インプレースおよびアウトオブプレースの転置/コピーマトリックスを示唆しています。 ドキュメントへのリンクです。 より速い10のインプレースで、最新バージョンのmklのドキュメントにいくつかの間違いが含まれているため、アウトオブプレースの実装を試すことをお勧めします。
配列のサイズが以前にわかっている場合は、組合を使用して支援することができます。 このような。
``。
namespace std;を使用します。
ユニオンua {。 int arr [2][3];。 int br [3][2];。 };。
int main(){。 ユニオンua uav;。 int karr [2][3] = {{{1,2,3}、{4,5,6}}};。 memcpy(uav.arr、karr、sizeof(karr));。 (int i = 0; i< 3; i ++)。 {。 (int j = 0; j< 2; j ++)。 cout<< uav.brr [i] [j]<< "";。 cout<< '\ n';。 }。
戻り0;。
}。 ``。
最新の線形代数ライブラリには、最も一般的な演算の最適化バージョンが含まれています。 それらの多くには、動的CPUディスパッチが含まれています。これは、プログラムの実行時にハードウェアに最適な実装を選択します(移植性を損なうことなく)。
これは通常、ベクトル拡張の固有の機能を介してファンクティーノを手動で最適化するためのより良い代替手段です。 後者は、実装を特定のハードウェアベンダーとモデルに結び付けます。別のベンダーにスワップする場合(例:. Power、ARM)または新しいベクトル拡張機能(例:. AVX512)、それらを最大限に活用するには、再度実装する必要があります。
たとえば、MKL転置には、BLAS拡張機能「imatcopy」が含まれます。 OpenBLASなどの他の実装でも見つけることができます。 ``。
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 ++:を利用できます。 ``。
void transpose(arma :: mat& matrix){。 arma :: inplace_trans(matrix);。 }。 ``。
私は、最も速い方法はO(n ^ 2)よりも高くてはいけないと思います。この方法でも、O(1)スペースだけを使用できます。 これを行う方法は、マトリックスを転置するときに実行することが次のとおりであるため、ペアで交換することです。M [i] [j] = M [j] [i]、M [i] [j]を一時的に格納し、次にM [i] [j] = M [j] [i]、および最後のステップ:M [j] [i] = temp。 これは1つのパスで行うことができるため、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;
}