[解決済み] 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

一番早い方法は何ですか？

どのように解決するのですか？

これは良い質問です。 例えば、行列の乗算やガウススミアなど、座標を交換するだけでなく、メモリ内で実際に行列を転置したいと思う多くの理由があります。

まず、私が転置のために使用する関数の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

最後に、私が行列の乗算で実際に行うことは（そしてガウス・スミアで）、正確に転置を取るのではなく、あるベクトルサイズ（例えば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

編集します。

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