difference result of mkl_sparse_s_export_csr in two subroutines

Tianxiong_Lu · ‎09-16-2015

Hi all,

I am using mkl_sparse_s_export_csr function for reading data in sparse_matrix_t object, but the return result is error while in other function. Here is the example:

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

void function1(const sparse_matrix_t csrA)
{
    // read data from sparse matrix
    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);
    printf("Matrix row count in function1: %d\n", row);
}

int main()
{
    //*******************************************************************************
    //     Declaration and initialization of parameters for sparse representation of
    //     the matrix A in the compressed sparse row format:
    //*******************************************************************************
#define M 5
#define N 5
#define NNZ 13
    //*******************************************************************************
    //    Sparse representation of the matrix A
    //*******************************************************************************
    float csrVal[NNZ]    = { 1.0, -1.0,     -3.0,
                             -2.0,  5.0,
                             4.0, 6.0, 4.0,
                             -4.0,       2.0, 7.0,
                             8.0,          -5.0 };
    MKL_INT    csrColInd[NNZ] = { 0,      1,        3,
                                  0,      1,
                                  2,   3,   4,
                                  0,           2,   3,
                                  1,             4 };
    MKL_INT    csrRowPtr[M+1] = { 0, 3, 5, 8, 11, 13 };
    // // Structure with sparse matrix stored in CSR format
    sparse_matrix_t       csrA;

    // Create handle with matrix stored in CSR format
    mkl_sparse_s_create_csr ( &csrA, SPARSE_INDEX_BASE_ZERO,
                              N,  // number of rows
                              M,  // number of cols
                              csrRowPtr,
                              csrRowPtr+1,
                              csrColInd,
                              csrVal );

    // read data from sparse matrix
    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);
    printf("Matrix row count in main: %d\n", row);

    // read data from function
    function1(csrA);

    mkl_sparse_destroy(csrA);
    return 1;
}

Result of this program is :

Matrix row count in main: 5
Matrix row count in function1: 0

Is there any error in function1? Why these two result is different?

thanks.

Tianxiong Lu

Ying_H_Intel · ‎09-16-2015

Hi TianXiong,

I can't reproduce the problem with both x64 and ia32 version. It can run fine in MSVC environments.

Best Regards,

Ying

Ying_H_Intel · ‎09-16-2015

Just for forum user's reference.

The root cause of the problem is the ILP64 library was linked.

Please note:

ILP64: 64bit integer

LP64: 32bit integer

MKL support both mode. But if one hope to use ILP64, please make sure all required integer was 64bit, for example, define as long long or use MKL_INT instead int type.

Best Regards,
Ying

mecej4 · ‎09-17-2015

Most of the MKL routines that deal with sparse matrices use the three-array CSR representation, but the routines that you have used here use the four-array CSR representation, which contains some redundant information. Furthermore, I find the description of the pointerB and pointerE arrays a bit confusing. For example, in the section on Sparse Matrix Storage Formats of the MKL reference manual we find "Element j of this integer array gives the index". This is ambiguous, especially to a C programmer. Do the element numbers j begin with 0 or 1? We have to look at the example in the manual to resolve this doubt. Next, take the description of the array pointerE: "An integer array that contains row indices, such that pointerE(j)-pointerB(1) is the index of the element in the values array that is last non-zero element in a row j of A." This is a needlessly complicated definition, and the names "pointerB" "pointerE" suggest that we may be talking about the Beginning and End of the elements in a particular row, but that is not the case, since pointerE gives us the index of the beginning of the next row, i.e., row j+1, instead of the end of row j.

I realize that this convention comes to us from NIST, but the confusion is real (at least for me!). Now that Intel has separated the MKL manual into C and Fortran versions, it would be nice to take the next step, and make all text concerning 2-D arrays follow the row-major and 0-based convention in the C version and, column-major and 1-based in the Fortran version.

In the specific example, Ying H has explained part of the reason. Here is what happens if the pointer array contains 8-byte integers when 4-byte integers are expected. The values [0, 3, 5, 8, 11, 13], stored as 8-byte integers will be seen as [0, 0, 3, 0, 5, 0] by the routine that expects 4-byte integers. That explains why one of the matrix dimensions is reported as 0.

I admit that it is possible that I am viewing this matter without my glasses on, and I should be grateful to have the correct interpretation pointed out. Part of the confusion could be caused by leaving out some parts of the NIST report http://www.netlib.org/utk/papers/sparse.ps when writing the MKL manual sections. The NIST report, for example, describes how the 4-array representation gives more flexibility in representing triangular and other matrices.

Ying_H_Intel · ‎09-17-2015

Hi Mecej4,

Thank you a lot for your feedback. I will forward them to our developer to see if there any improvements.

Your observation are exact right, There are 3 array CSR representation and 4- representation. MKL use them both.

The 3-array CSR was often used in Direct sparse solvers, PARSDIO and DSS, and some of sparse matrix format conversion.

The 4-array CSR is from NIST Sparse BLAS library, been used in Sparse BLAS Levels 2 and Level 3, more often in latest Inspector-executor Sparse BLAS (SpMV 2) Routines. for example, this function

The 3-array variation of the CSR format has a restriction: all non-zero elements are stored continuously, that is the set of non-zero elements in the row J goes just after the set of non-zero elements in the row J-1.

There are no such restrictions in the general (NIST) CSR format. This may be useful, for example, if there is a need to operate with different submatrices of the matrix at the same time. In this case, it is enough to define the arrays pointerB and pointerEfor each needed submatrix so that all these arrays are pointers to the same array values.

Regarding that row-major and 0-based convention in the C version and, column-major and 1-based in the Fortran version. right , it should valid for most of MKL fucntion, But the exception is the Sparse matrix, like here SPARSE BLAS 2 and 3, new SPMV 2 interface, as i recalled, we have get cross requests include 0-based, column-major, 1 based, row-major etc. both are supported in these API.

Best Regards,

Ying

Tianxiong_Lu · ‎09-18-2015

Now I am confused by the 4-array CSR representation in MKL. I get these informations in MKL 11.3 reference manual and example code. Which is wrong and which is correct?

1. In description of routine mkl_sparse_?_create_csr, https://software.intel.com/zh-cn/node/590110

rows_end

Array of at least length m. This array contains row indices, such that rows_end[i] - rows_start[0] - 1 is the last index of row i in the arrays values and col_indx.

Refer to pointerE array description in CSR Format for more details.

It means that rows_end is the first index of next row.

2. In description of Sparse Matrix Storage Formats, https://software.intel.com/zh-cn/node/522243#MKL_APPA_SMSF_3

pointerE: An integer array that contains row indices, such that pointerE[j]-pointerB[0] is the index of the element in the values array that is last non-zero element in a row j of A

In this page, it means that pointerE is the last index of this row.

3. In the example file of MKL , in directory spblasc. The rows_end indices is case 1, for example:

//*******************************************************************************

    //    Sparse representation of the matrix A
    //*******************************************************************************
    double csrVal[NNZ]    = { 1.0, -1.0,     -3.0,
                             -2.0,  5.0,
                                         4.0, 6.0, 4.0,
                             -4.0,       2.0, 7.0,
                                    8.0,          -5.0 };
    MKL_INT    csrColInd[NNZ] = { 0,      1,        3,
                              0,      1,
                                           2,   3,   4,
                              0,           2,   3,
                                      1,             4 };
    MKL_INT    csrRowPtr[M+1] = { 0, 3, 5, 8, 11, 13 };
    // Descriptor of main sparse matrix properties
    struct matrix_descr descrA;
    // // Structure with sparse matrix stored in CSR format
    sparse_matrix_t       csrA;
    //*******************************************************************************
    //    Declaration of local variables:
    //*******************************************************************************
    double x_m[M*NRHS]  = { 1.0, 5.0, 3.0, 4.0, 2.0, 2.0, 10.0, 6.0, 8.0, 4.0};
    double y_m[M*NRHS]  = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    double x_v  = { 3.0, 2.0, 5.0, 4.0, 1.0};
    double y_v  = { 0.0, 0.0, 0.0, 0.0, 0.0};
    double tmp_v  = { 0.0, 0.0, 0.0, 0.0, 0.0};
    double alpha = 1.0, beta = 0.0;
    MKL_INT    i;

    printf( "\n EXAMPLE PROGRAM FOR CSR format routines from package\n" );
    printf( "-------------------------------------------------------\n" );

    // Create handle with matrix stored in CSR format
    mkl_sparse_d_create_csr ( &csrA, SPARSE_INDEX_BASE_ZERO,
                                    M,  // number of rows
                                    M,  // number of cols
                                    csrRowPtr,
                                    csrRowPtr+1,
                                    csrColInd,
                                    csrVal );