4. 稀疏矩阵转置的 CSR 实现与性能分析

CSR 转置的核心结构与两阶段代价

CSR 用 row_ptr（长度 m+1）划分行、用 col_idx 与 vals 按行内顺序存放非零元。转置后得到的是 Aᵀ 的 CSR，等价于把 A 的每个非零元 (i, j, v) 变成 (j, i, v) 并按新行（原列）分组。直接逐元素插入会破坏行内顺序且需要动态扩容，因此标准做法是计数排序：先统计 A 每一列的非零个数得到 Aᵀ 的行长度，前缀和得到 Aᵀ 的 row_ptr，再用一个游标数组把每个元素放到 Aᵀ 对应行的下一个空位。

参考实现（C++，CSR 三元组结构）

下面给出可直接编译运行的完整实现，包含转置函数与对称性验证。矩阵采用 0 基索引，row_ptr 长度为 m+1。

#include <cstdio>
#include <vector>
#include <cassert>
using namespace std;

struct CSR {
    int m, n, nnz;
    vector<int> row_ptr;   // 长度 m+1
    vector<int> col_idx;   // 长度 nnz
    vector<double> vals;   // 长度 nnz
};

// 计数排序实现 CSR 转置：返回 A^T 的 CSR
CSR transpose_csr(const CSR& A) {
    CSR T;
    T.m = A.n; T.n = A.m; T.nnz = A.nnz;
    T.row_ptr.assign(T.m + 1, 0);
    T.col_idx.assign(T.nnz, 0);
    T.vals.assign(T.nnz, 0.0);

    // 阶段 1：按列计数
    for (int k = 0; k < A.nnz; ++k)
        T.row_ptr[A.col_idx[k] + 1]++;
    // 前缀和
    for (int j = 0; j < T.m; ++j)
        T.row_ptr[j + 1] += T.row_ptr[j];

    // 阶段 2：游标填充
    vector<int> next(T.row_ptr.begin(), T.row_ptr.end() - 1);
    for (int i = 0; i < A.m; ++i) {
        for (int k = A.row_ptr[i]; k < A.row_ptr[i + 1]; ++k) {
            int j = A.col_idx[k];
            int pos = next[j]++;
            T.col_idx[pos] = i;
            T.vals[pos] = A.vals[k];
        }
    }
    return T;
}

// 稠密展开用于对照
vector<vector<double>> to_dense(const CSR& A) {
    vector<vector<double>> D(A.m, vector<double>(A.n, 0.0));
    for (int i = 0; i < A.m; ++i)
        for (int k = A.row_ptr[i]; k < A.row_ptr[i + 1]; ++k)
            D[i][A.col_idx[k]] = A.vals[k];
    return D;
}

int main() {
    // 3x4 示例矩阵
    CSR A;
    A.m = 3; A.n = 4; A.nnz = 6;
    A.row_ptr = {0, 2, 4, 6};
    A.col_idx = {0, 2, 1, 3, 0, 2};
    A.vals    = {5, 7, 6, 8, 9, 4};

    CSR T = transpose_csr(A);

    // 结构检查
    assert(T.m == 4 && T.n == 3 && T.nnz == 6);
    assert(T.row_ptr.back() == A.nnz);

    // 稠密对照：T 应等于 A 的转置
    auto DA = to_dense(A);
    auto DT = to_dense(T);
    for (int i = 0; i < A.m; ++i)
        for (int j = 0; j < A.n; ++j)
            assert(DA[i][j] == DT[j][i]);

    // 对称性：转置两次还原
    CSR TT = transpose_csr(T);
    assert(TT.m == A.m && TT.n == A.n && TT.nnz == A.nnz);
    for (int i = 0; i < A.m; ++i)
        for (int k = A.row_ptr[i]; k < A.row_ptr[i + 1]; ++k) {
            int j = A.col_idx[k];
            bool found = false;
            for (int q = TT.row_ptr[i]; q < TT.row_ptr[i + 1]; ++q)
                if (TT.col_idx[q] == j && TT.vals[q] == A.vals[k]) found = true;
            assert(found);
        }

    printf("转置与对称性验证通过\n");
    return 0;
}

复杂度分析与性能要点

关键点：填充阶段对 Aᵀ 的写入位置由 next[j] 决定，属于按列分散的随机写，缓存局部性弱于计数阶段。若矩阵列数 n 很大而每列非零数很少，随机写会显著拖慢性能；此时可考虑分块转置（按列块切分，使每块的目标行段落在较小范围内）以改善缓存命中。计数阶段只读 col_idx 并写计数数组，通常更快。