C语言实现数组转置的代码详解

  地址:    ... | +0     | +1     | ... | +n-1   | +n     | +n+1  | ...
  元素:    ... | a[0][0]| a[0][1]| ... | a[0][n-1]| a[1][0]| a[1][1]| ...

for (i = 0; i < n; ++i)
  for (j = i+1; j < n; ++j)
    swap(a[i*n + j], a[j*n + i]);

for (ii = 0; ii < rows; ii += b)
  for (jj = 0; jj < cols; jj += b)
    // 对矩阵子块 (ii..ii+b-1, jj..jj+b-1) 进行单独转置
    for (i = ii; i < min(ii+b, rows); ++i)
      for (j = jj; j < min(jj+b, cols); ++j)
        dst[j*rows + i] = src[i*cols + j];

#pragma omp parallel for collapse(2)
for (i = 0; i < rows; ++i)
  for (j = 0; j < cols; ++j)
    dst[j*rows + i] = src[i*cols + j];

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <omp.h>
 
#ifndef min
#define min(a,b) (((a)<(b))?(a):(b))
#endif
 
/**
 * 基本额外空间转置
 * @param src   原矩阵指针
 * @param rows  行数
 * @param cols  列数
 * @param dst   目标矩阵指针（已分配 rows*cols 大小）
 */
void transpose_with_buffer(int *src, int rows, int cols, int *dst) {
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            dst[j * rows + i] = src[i * cols + j];
        }
    }
}
 
/**
 * 方阵原地转置
 * 适用于 n x n 方阵
 */
void transpose_inplace(int *a, int n) {
    for (int i = 0; i < n; ++i) {
        for (int j = i + 1; j < n; ++j) {
            int tmp = a[i * n + j];
            a[i * n + j] = a[j * n + i];
            a[j * n + i] = tmp;
        }
    }
}
 
/**
 * 块转置 (块大小 block_size)
 * @param src         原矩阵
 * @param rows,cols   原矩阵尺寸
 * @param block_size  块大小
 * @param dst         目标矩阵
 */
void transpose_block(int *src, int rows, int cols, int block_size, int *dst) {
    for (int ii = 0; ii < rows; ii += block_size) {
        for (int jj = 0; jj < cols; jj += block_size) {
            int max_i = min(ii + block_size, rows);
            int max_j = min(jj + block_size, cols);
            for (int i = ii; i < max_i; ++i) {
                for (int j = jj; j < max_j; ++j) {
                    dst[j * rows + i] = src[i * cols + j];
                }
            }
        }
    }
}
 
/**
 * openmp 并行转置 (基本方法)
 */
void transpose_omp(int *src, int rows, int cols, int *dst) {
    #pragma omp parallel for collapse(2)
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            dst[j * rows + i] = src[i * cols + j];
        }
    }
}
 
/**
 * 性能测试主函数
 */
int main(int argc, char *argv[]) {
    int rows = 4096, cols = 4096;
    int *a = (int*)malloc(sizeof(int) * rows * cols);
    int *b = (int*)malloc(sizeof(int) * rows * cols);
    if (!a || !b) {
        fprintf(stderr, "内存分配失败\n");
        return exit_failure;
    }
 
    // 初始化
    for (int i = 0; i < rows; ++i)
        for (int j = 0; j < cols; ++j)
            a[i * cols + j] = i * cols + j;
 
    struct timespec t1, t2;
    double elapsed;
 
    // 1. 额外空间转置
    clock_gettime(clock_monotonic, &t1);
    transpose_with_buffer(a, rows, cols, b);
    clock_gettime(clock_monotonic, &t2);
    elapsed = (t2.tv_sec - t1.tv_sec) + (t2.tv_nsec - t1.tv_nsec)/1e9;
    printf("buffer transpose: %.6f s\n", elapsed);
 
    // 2. 原地方阵转置 (只针对方阵 a)
    clock_gettime(clock_monotonic, &t1);
    transpose_inplace(a, cols);
    clock_gettime(clock_monotonic, &t2);
    elapsed = (t2.tv_sec - t1.tv_sec) + (t2.tv_nsec - t1.tv_nsec)/1e9;
    printf("in-place square transpose: %.6f s\n", elapsed);
 
    // 3. 块转置
    clock_gettime(clock_monotonic, &t1);
    transpose_block(a, rows, cols, 64, b);
    clock_gettime(clock_monotonic, &t2);
    elapsed = (t2.tv_sec - t1.tv_sec) + (t2.tv_nsec - t1.tv_nsec)/1e9;
    printf("block transpose (64): %.6f s\n", elapsed);
 
    // 4. openmp 并行转置
    clock_gettime(clock_monotonic, &t1);
    transpose_omp(a, rows, cols, b);
    clock_gettime(clock_monotonic, &t2);
    elapsed = (t2.tv_sec - t1.tv_sec) + (t2.tv_nsec - t1.tv_nsec)/1e9;
    printf("openmp parallel transpose: %.6f s\n", elapsed);
 
    free(a);
    free(b);
    return 0;
}