mkl_sparse_sp2m Issue C++

[streamsky]

I'm trying to perform multiplication between two sparse matrices using mkl_sparse_sp2m and I'm not getting the expected result.  Below is my example code showing the computed and expected result. What am I doing wrong?:

#include <iostream>
#include <assert.h> 
#include "mkl_spblas.h"

void checkStat(const int& stat) {
    if(stat != SPARSE_STATUS_SUCCESS) {
        std::cout << "Status: " << stat << std::endl;
        assert(false);
    }
}

void outputCsr(sparse_matrix_t& A) {
    MKL_INT rowsC, colsC;
    sparse_index_base_t indexing = SPARSE_INDEX_BASE_ZERO;
    MKL_INT *cols_start = NULL, *cols_end = NULL, *row_indx = NULL;
    double *values_C = NULL;
    checkStat(mkl_sparse_d_export_csr(A,
        &indexing,
        &rowsC,
        &colsC,
        &cols_start,
        &cols_end,
        &row_indx,
        &values_C));

    int len_row_index = cols_end[colsC - 1] - cols_start[0];
    int k = rowsC;

    std::cout << "Rows, Cols = " << rowsC << ", " << colsC << std::endl;

    for(int i = 0; i < len_row_index; ++i)
        std::cout << values_C[i] << ", ";
    std::cout << std::endl;

    for(int i = 0; i < len_row_index; ++i)
        std::cout << row_indx[i] << ", ";
    std::cout << std::endl;

    for(int i = 0; i < k; ++i)
        std::cout << cols_start[i] << ", ";
    std::cout << std::endl;

    for(int i = 0; i < k; ++i)
        std::cout << cols_end[i] << ", ";
    std::cout << std::endl;
}

int main() {
    int m = 3;
    int k = 3;

    sparse_matrix_t A;

    /* Matrix in CRS format
    *
    * { { 0,  0,  1 } 
    *   { 0, -1,  2 } 
    *   { 1,  0,  0 } } 
    */
    double vals[] = { 1, -1, 2, 1 };
    int cols[]   = { 2, 1, 2, 0 };
    int ptrB[]  = { 0, 1, 3 };
    int ptrE[]  = { 1, 3, 4 };

    checkStat(mkl_sparse_d_create_csr(&A,
        SPARSE_INDEX_BASE_ZERO, // indexing (base),
        m,      // rows,
        k,      // cols,
        ptrB,    // rows_start,
        ptrE,    // rows_end,
        cols,      // col_indx,
        vals));
    std::cout << "A" << std::endl;
    outputCsr(A);

    matrix_descr descr;
    descr.type = SPARSE_MATRIX_TYPE_GENERAL;

    sparse_matrix_t B;

    checkStat(mkl_sparse_sp2m(
        SPARSE_OPERATION_NON_TRANSPOSE,
        descr,
        A,
        SPARSE_OPERATION_NON_TRANSPOSE,
        descr,
        A,
        SPARSE_STAGE_FULL_MULT,
        &B));

    std::cout << "Result:" << std::endl;
    outputCsr(B);

    /* Expected answer
    *
    * { { 1,  0,  0 } 
    *   { 2,  1,  -2 } 
    *   { 0,  0,  1 } } 
    * 
    *  vals: 1, 2, 1, -2
    *  cols: 0, 0, 1, 2, 2
    *  ptrB: 0, 1, 4,
    *  prtE: 1, 2, 5
    * 
    * MKL Answer
    * { {1, 0, 0}
    *   {0, 1, -2}
    *   {2, 0, 1}}
    * 
    *   vals: 1, 1, -2, 2, 1, 
    *   cols: 0, 1, 2, 0, 2, 
    *   ptrB: 0, 1, 4, 
    *   ptrE: 1, 4, 5
    */

    checkStat(mkl_sparse_destroy(A));

    return 0;
}

[VidyalathaB_Intel · Accepted solution]

Hi David,

Thanks for reaching out to us.

The output column/row indices for sp2m are not sorted, this is implied by the algorithm currently used.

In order to sort the indices, you can call mkl_sparse_order for the resultant matrix (here it is B) after using mkl_sparse_sp2m.

Just insert the mkl_sparse_order function as below in your code 

std::cout << "Result:" << std::endl;
checkStat(mkl_sparse_order (B));
outputCsr(B);

After making the above changes, recompile the code and you will get the expected results.

Please do let us know if you face any issues.

Regards,

Vidya.

[streamsky]

Vidya,

Thank you for the solution.   Could you please expand upon "implied by the algorithm" and the use of mkl_sparse_order?  I haven't found these insights in the developer reference.  For example, does mkl_sparse_order need to be called to view a matrix in CSR format after MKL sparse matrix operations?

Also, for anyone else reading this, the expect answer should have been:

vals: 1, 2, 1, -2, 1

cols: 0, 0, 1, 2, 2

ptrB: 0, 1, 4,

prtE: 1, 4, 5

David

[Kirill_V_Intel]

Hi David,

We'll update the documentation of sp2m to say about potential unsortedness.

If you are curious about the reasons: while oneMKL neither exposes any information about the algorithm used or makes any promises to use only one of the existing ones, a very well known Gustavson algorithm for sparse * sparse matrix multiplication naturally produces unsorted indices for the output.

As a lot of functionality actually does not need the indices to be sorted, it would be an unnecessary performance penalty for these cases to do the sorting additionally within sp2m (by default, unless a special knob for sorted/unsorted output appears).

Best,
Kirill

[VidyalathaB_Intel]

Hi David,

Thanks for accepting our solution.

As your issue is resolved we are going ahead and closing this thread. Please post a new question if you require any additional assistance from Intel as this thread will no longer be monitored.

Have a Nice Day!

Regards,

Vidya.