A strange result of sparse matrix addition with mol_sparse_s_add

[Tianxiong_Lu]

Hi all，

I have a new question about sparse matrix addition routine mkl_sparse_s_add, it return a strange result. And when I using double precision routine ml_sparse_d_add, everything is ok. 

Following is the example program.

#include "mkl_spblas.h"
#include <stdio.h>

void print_csr_sparse_s(const sparse_matrix_t csrA)
{
    // Read Matrix Data and Print it
    int row, col;
    sparse_index_base_t indextype;
    int * bi, *ei;
    int * j;
    float *rv;
    sparse_status_t status = mkl_sparse_s_export_csr(csrA, &indextype, &row, &col, &bi, &ei, &j, &rv);
    if (status==SPARSE_STATUS_SUCCESS)
    {
        printf("SparseMatrix(%d x %d) [base:%d]\n", row, col, indextype);
        for (int r = 0; r<row; ++r)
        {
            for (int idx = bi; idx<ei; ++idx)
            {
                printf("  <%d, %d> \t %f\n", r, j[idx], rv[idx]);
            }
        }
    }
    return;
}

void print_csr_sparse_d(const sparse_matrix_t csrA)
{
    // Read Matrix Data and Print it
    int row, col;
    sparse_index_base_t indextype;
    int * bi, *ei;
    int * j;
    double *rv;
    sparse_status_t status = mkl_sparse_d_export_csr(csrA, &indextype, &row, &col, &bi, &ei, &j, &rv);
    if (status==SPARSE_STATUS_SUCCESS)
    {
        printf("SparseMatrix(%d x %d) [base:%d]\n", row, col, indextype);
        for (int r = 0; r<row; ++r)
        {
            for (int idx = bi; idx<ei; ++idx)
            {
                printf("  <%d, %d> \t %f\n", r, j[idx], rv[idx]);
            }
        }
    }
    return;
}

// test addition of sparse matrix
void test_add_s()
{
    // Define sparse-matrix M
    int mi[5] = {0, 2, 5, 8, 10};
    int mj[10] = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3};
    float mv[10] = {2.0f, 1.0f, 1.0f, 2.0f, 1.0f, 1.0f, 2.0f, 1.0f, 1.0f, 2.0f};
    sparse_matrix_t M;

    // Define sparse-matrix N
    int ni[5] = {0, 1, 2, 3, 4};
    int nj[4] = {0, 1, 2, 3};
    float nv[4] = {3.0f, 2.0f, 1.0f, -1.0f};
    sparse_matrix_t N;

    // create csr matrix
    mkl_sparse_s_create_csr(&M, SPARSE_INDEX_BASE_ZERO, 4, 4, mi, mi+1, mj, mv);
    mkl_sparse_s_create_csr(&N, SPARSE_INDEX_BASE_ZERO, 4, 4, ni, ni+1, nj, nv);
    // output matrix
    print_csr_sparse_s(M);
    print_csr_sparse_s(N);

    // do addition
    sparse_matrix_t C;
    mkl_sparse_s_add(SPARSE_OPERATION_NON_TRANSPOSE ,M, 2, N, &C);

    // output result
    print_csr_sparse_s(C);

    // free memory
    mkl_sparse_destroy(M);
    mkl_sparse_destroy(N);
    mkl_sparse_destroy(C);
}

void test_add_d()
{
    // Define sparse-matrix M
    int mi[5] = {0, 2, 5, 8, 10};
    int mj[10] = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3};
    double mv[10] = {2.0f, 1.0f, 1.0f, 2.0f, 1.0f, 1.0f, 2.0f, 1.0f, 1.0f, 2.0f};
    sparse_matrix_t M;

    // Define sparse-matrix N
    int ni[5] = {0, 1, 2, 3, 4};
    int nj[4] = {0, 1, 2, 3};
    double nv[4] = {3.0f, 2.0f, 1.0f, -1.0f};
    sparse_matrix_t N;

    // create csr matrix
    mkl_sparse_d_create_csr(&M, SPARSE_INDEX_BASE_ZERO, 4, 4, mi, mi+1, mj, mv);
    mkl_sparse_d_create_csr(&N, SPARSE_INDEX_BASE_ZERO, 4, 4, ni, ni+1, nj, nv);
    // output matrix
    print_csr_sparse_d(M);
    print_csr_sparse_d(N);

    // do addition
    sparse_matrix_t C;
    mkl_sparse_d_add(SPARSE_OPERATION_NON_TRANSPOSE ,M, 2, N, &C);

    // output result
    print_csr_sparse_d(C);

    // free memory
    mkl_sparse_destroy(M);
    mkl_sparse_destroy(N);
    mkl_sparse_destroy(C);
}

int main()
{
    test_add_d();
	test_add_s();
    return 1;
}

I got results twice, the first is correct, and the second is wrong.

SparseMatrix(4 x 4) [base:0]
  <0, 0>         2.000000
  <0, 1>         1.000000
  <1, 0>         1.000000
  <1, 1>         2.000000
  <1, 2>         1.000000
  <2, 1>         1.000000
  <2, 2>         2.000000
  <2, 3>         1.000000
  <3, 2>         1.000000
  <3, 3>         2.000000
SparseMatrix(4 x 4) [base:0]
  <0, 0>         3.000000
  <1, 1>         2.000000
  <2, 2>         1.000000
  <3, 3>         -1.000000
SparseMatrix(4 x 4) [base:0]
  <0, 0>         7.000000
  <0, 1>         2.000000
  <1, 0>         2.000000
  <1, 1>         6.000000
  <1, 2>         2.000000
  <2, 1>         2.000000
  <2, 2>         5.000000
  <2, 3>         2.000000
  <3, 2>         2.000000
  <3, 3>         3.000000
SparseMatrix(4 x 4) [base:0]
  <0, 0>         2.000000
  <0, 1>         1.000000
  <1, 0>         1.000000
  <1, 1>         2.000000
  <1, 2>         1.000000
  <2, 1>         1.000000
  <2, 2>         2.000000
  <2, 3>         1.000000
  <3, 2>         1.000000
  <3, 3>         2.000000
SparseMatrix(4 x 4) [base:0]
  <0, 0>         3.000000
  <1, 1>         2.000000
  <2, 2>         1.000000
  <3, 3>         -1.000000
SparseMatrix(4 x 4) [base:0]
  <0, 0>         3.000000
  <0, 1>         0.000000
  <1, 0>         0.000000
  <1, 1>         2.000000
  <1, 2>         0.000000
  <2, 1>         0.000000
  <2, 2>         1.000000
  <2, 3>         0.000000
  <3, 2>         0.000000
  <3, 3>         -1.000000

Best regrads,

Tianxiong Lu

[mecej4]

You have the same sort of errors as in your other thread, please see https://software.intel.com/en-us/comment/reply/593560/1839991?quote=1#comment-form . In the 4-array version of the CSR sparse matrix representation, you have to supply the indices of the beginning of each row and the end of that row. In the 3-array version, row i is assumed to contain items ia to ia[i+1]-1. In the 4-array version, row i is assumed to contain items ia_start to ia_end-1.

[Tianxiong_Lu]

HI mecej4,

I don't think that the end_row indices array is error, because the double precision version is correct.

From single precision version's result, we can know the sparse matrix structure is created for addition, but the result only include part of B, alpha * op(A) is missing.

best regards

Tianxiong Lu

[Ying_H_Intel]

Hi tianxiong, 

The indices array should be fine.  I can reproduce the problem. It seems a bug in single add function. The issue have been forwarded to the deveoper. We uioo update you if any updates. 

Best Regards,
Ying