这是个好问题。有很多原因您希望在内存中实际转置矩阵，而不仅仅是交换坐标，例如在矩阵乘法和高斯拖尾。
首先让我列出我用于转置的一个函数（编辑：请参阅我的答案末尾，在那里我找到了一个更快的解决方案）
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
这是英特尔的一篇论文，解释说
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];
}
}
编辑：
我尝试了几个函数来为大型矩阵找到最快的转置。最后，最快的结果是使用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);
}
}
}
}
}