Search code examples
c++matlabcudaopenmpthrust

Find the row sum of a matrix using CSR or COO storing format


I have for example the following matrix B which is stored in COO and CSR format (retrieved from the non-symetric example here). Could you please suggest an efficient c++ way to apply the matlab sum(B,2) function using the coo or csr(or both) storing format? Because it is quit possible to work with large arrays can we do that using parallel programming (omp or CUDA (e.g, thrust))?

Any algorithmic or library based suggestions are highly appreciated. Thank you!

PS: Code to construct a sparse matrix and get the CSR coordinates can be found for example in the answer of this post.

enter image description here

COO format: CSR format:

       row_index col_index value                 columns  row_index  value
          1         1         1                     0         0       1
          1         2        -1                     1         3      -1
          1         3        -3                     3         5      -3   
          2         1        -2                     0         8      -2
          2         2         5                     1         11      5
          3         3         4                     2         13      4
          3         4         6                     3                 6
          3         5         4                     4                 4
          4         1        -4                     0                -4
          4         3         2                     2                 2
          4         4         7                     3                 7     
          5         2         8                     1                 8
          5         5        -5                     4                -5

Solution

  • For COO its pretty simple:

    struct MatrixEntry {
        size_t row;
        size_t col;
        int value;
    };
    
    std::vector<MatrixEntry> matrix = {
        { 1, 1, 1 },
        { 1, 2, -1 },
        { 1, 3, -3 },
        { 2, 1, -2 },
        { 2, 2, 5 },
        { 3, 3, 4 },
        { 3, 4, 6 },
        { 3, 5, 4 },
        { 4, 1, -4 },
        { 4, 3, 2 },
        { 4, 4, 7 },
        { 5, 2, 8 },
        { 5, 5, -5 },
    };
    
    std::vector<int> sum(5);
    
    for (const auto& e : matrix) {
        sum[e.row-1] += e.value;
    }
    

    and for large matrixes you can just split up the for loop into multiple smaller ranges and add the results at the end.

    If you only need the sum of each row (and not columwise) CSR is also straight forward (and even more efficient):

    std::vector<int> row_idx = { 0, 3, 5, 8, 11, 13 };
    std::vector<int> value = { 1, -1, -3, -2, 5, 4, 6, 4, -4, 2, 7, 8, -5 };
    
    std::vector<int> sum(5);
    
    for(size_t i = 0; i < row_idx.size()-1; ++i) {
        sum[i] = std::accumulate(value.begin() + row_idx[i], value.begin() + row_idx[i + 1], 0);
    }
    

    Again, for parallelism you can simply split up the loop.