matrixcompressionrowsparse-matrixtranspose

Compressed Sparse Row Transpose


As you know we can write sparse matrices in compressed row storage (CRS) (or alternatively, compressed sparse row (CSR)). Let A be an m n matrix. The transpose of A is an n x m matrix A' such that for all 0 <= i < n and 0 <= j < m, A'(i; j) = A(j; i).

I need to write the algorithm for transposing a matrix in CRS representation. How can i approach this problem?


Solution

  • I was looking for something like that. Here is my algorithm. I don't know if it is the fastest, but I think it's quite good.

    EDIT: Essentially the same algorithm is implemented in C++ module for scipy.

    Suppose matrix is represent by this struct:

    struct CRSMatrix
    {
        int n; // number of rows
        int m; // number of columns
        int nz; // number of non-zero elements
        std::vector<double> val; // non-zero elements
        std::vector<int> colIndex; // column indices
        std::vector<int> rowPtr; // row ptr
    };
    

    This function does it:

    CRSMatrix sparse_transpose(const CRSMatrix& input) {
        CRSMatrix res{
            input.m,
            input.n,
            input.nz,
            std::vector<double>(input.nz, 0.0),
            std::vector<int>(input.nz, 0),
            std::vector<int>(input.m + 2, 0) // one extra
        };
    
        // count per column
        for (int i = 0; i < input.nz; ++i) {
            ++res.rowPtr[input.colIndex[i] + 2];
        }
    
        // from count per column generate new rowPtr (but shifted)
        for (int i = 2; i < res.rowPtr.size(); ++i) {
            // create incremental sum
            res.rowPtr[i] += res.rowPtr[i - 1];
        }
    
        // perform the main part
        for (int i = 0; i < input.n; ++i) {
            for (int j = input.rowPtr[i]; j < input.rowPtr[i + 1]; ++j) {
                // calculate index to transposed matrix at which we should place current element, and at the same time build final rowPtr
                const int new_index = res.rowPtr[input.colIndex[j] + 1]++;
                res.val[new_index] = input.val[j];
                res.colIndex[new_index] = i;
            }
        }
        res.rowPtr.pop_back(); // pop that one extra
    
        return res;
    }