CUDA N×N矩阵原地转置（仅遍历上三角）线程索引分配问题

我正在实现CUDA环境下N×N矩阵的原地转置，核心思路是仅遍历上三角元素（对角线元素转置后保持不变），通过单个线程用临时变量交换A[i,j]与A[j,i]的值。但当前代码中分配下一组待处理索引对A[i_next,j_next]的两行代码存在错误，无法让每个线程仅访问正确的矩阵元素。以N=5、总线程数=3的矩阵为例，需要修正这两行代码的逻辑。

原错误代码 #

__global__ void in_place_transpose(double *A, int N)
{
    int t_idx = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned int total_threads = blockDim.x * gridDim.x;
    int row_idx = t_idx / N;
    int col_idx = t_idx % N;

    // only access upper triangular entries
    if (row_idx < col_idx)
    {
        while (row_idx < (N - 1) && col_idx < N) // A[N-1,N] is the last entry we change
        {
            int to = row_idx * N + col_idx;
            int from = col_idx * N + row_idx;

            if (from != to) // diagonal is constant
            {
                // printf("from: %d, to: %d \n", from, to); // show operation pair
                // printf("th: %d, row: %d, col: %d \n", t_idx, row_idx, col_idx); // show indizes
                int temp = A[to];
                A[to] = A[from];
                A[from] = temp;
            }
            // assign next matrix entry to thread (FAULTY)
            row_idx = (int)((row_idx + total_threads) % (N - row_idx - 1));
            col_idx = (int)(col_idx + (total_threads % (N - row_idx - 1)) % (N - row_idx - 1));
            
        }
    }
}

问题分析 #

原索引更新逻辑完全错误，没有遵循上三角元素的遍历规则。正确的做法应该是：将所有上三角元素（row < col）按行优先顺序映射为线性序列，每个线程从初始分配的位置开始，每次跳过total_threads个元素，保证每个线程处理的元素不重复、不越界，且始终处于上三角区域。

修正后的代码 #

__global__ void in_place_transpose(double *A, int N)
{
    int t_idx = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned int total_threads = blockDim.x * gridDim.x;
    // 计算上三角元素的总数
    unsigned int upper_count = N * (N - 1) / 2;

    // 每个线程处理线性偏移量为 t_idx, t_idx+total_threads, ... 的上三角元素
    for (unsigned int k = t_idx; k < upper_count; k += total_threads)
    {
        // 将线性偏移量k转换为对应的(row, col)
        int row = 0;
        // 找到最大的row满足 row*(2*N - row -1)/2 <= k
        while ((row + 1) * (2 * N - (row + 1) - 1) / 2 <= k)
        {
            row++;
        }
        int col = row + 1 + (k - row * (2 * N - row - 1) / 2);

        int to = row * N + col;
        int from = col * N + row;

        // 交换元素（对角线元素无需处理，这里已经保证row<col，所以from!=to）
        double temp = A[to];
        A[to] = A[from];
        A[from] = temp;
    }
}

修正逻辑说明 #

1. 线性化上三角元素：上三角区域共有N*(N-1)/2个元素，按行优先顺序（先遍历row=0的所有col>0，再row=1的col>1，以此类推）映射为从0到upper_count-1的线性偏移量。
2. 线程分配规则：每个线程处理偏移量为t_idx, t_idx+total_threads, t_idx+2*total_threads...的元素，直到超出上三角元素总数。
3. 线性偏移转行列索引：通过数学计算找到偏移量k对应的row和col，确保始终满足row < col，避免访问下三角或对角线元素。

这种方式能保证每个线程只处理属于自己的上三角元素，不会重复或越界，完美适配原地转置的需求。

内容的提问来源于stack exchange，提问作者grube