c - 编写函数以使用英特尔 MKL 库进行稀疏矩阵乘法

Question

我正在尝试开发一个将矩阵 A 和 B 相乘的函数，它们是通用格式但本质上是稀疏的。这些矩阵包含复数。我的问题是，当我不使用该函数并将所有内容写入 main() 中时，乘法对于任何大小的数组都完美有效。但是当我使用自己的函数时，结果会损坏，并且大多数时候我会收到随机错误消息。

此函数的作用如下:

1. 将 A 和 B 转换为 CSR 格式 (mkl_zdnscsr)。
2. 使用第 1 步 (mkl_sparse_z_create_csr) 中的数据为 A 和 B 创建 CSR 句柄。
3. 使用步骤 2 中的句柄 (mkl_sparse_spmm) 将 A 和 B 相乘，并将输出存储在“结果”中。

我的猜测是我从函数返回“结果”的方式在某种程度上是不正确的，因为我已经检查了步骤 1 和 2 的输出并且它们产生了正确的输出。

知道问题是什么吗？我将在下面包含我的代码的摘要版本以供您引用。

提前非常感谢您。

阿夫辛

/* ***************** Macro ********************* */
#define ALIGN 128

/* To avoid constantly repeating the part of code that checks different functions' status, using the below macros */
#define CHECK_SPARSE(function) do { \
if(function != SPARSE_STATUS_SUCCESS)             \
 {                                                 \
 status = 2;                                       \
 goto memory_free;                                 \
 }                                                 \
} while(0)


/* ****************** Main ******************** */
int main()
{

 << matrices A and B are generated using some data >>

    MKL_INT stat = 0;

    // This part calls the function to multiply matrices as I discussed.

    // A x B --> csrC
    sparse_matrix_t  csrC = NULL;
    stat = dnmm_sp_CSR_handle(CfPrime, Num_of_Buses, Num_of_Branches, CfPrime_nonzero, Yf, Num_of_Branches, Num_of_Buses, Yf_nonzero, &csrC);
    printf("\nstat = %i", stat);

    // Now I convert the csrC to 4-array version of CSR.
    MKL_INT rows, cols;
    sparse_index_base_t indexing = 0;
    MKL_INT *columns_C = NULL, *pointerB_C = NULL, *pointerE_C = NULL;
    MKL_Complex16  *values_C = NULL;
    mkl_sparse_z_export_csr(csrC, &indexing, &rows, &cols, &pointerB_C, &pointerE_C, &columns_C, &values_C);

    // Print the number of rows and columns of converted matrix (which are incorrect sizes)
    printf("\nrows = %i , cols = %i", rows, cols);
}

/* ****************** Function ******************** */
// This function receives two dense matrice, convert them to sparse CSR format, multiply them, and returns the result in CSR handle
int dnmm_sp_CSR_handle(MKL_Complex16 *A, MKL_INT A_rownum, MKL_INT A_colnum, MKL_INT A_nnz, MKL_Complex16 *B, MKL_INT B_rownum, MKL_INT B_colnum, MKL_INT B_nnz, sparse_matrix_t *result) {

    // A : Matrix A
    // A_rownum : Number of rows in matrix A
    // A_colnum : Number of columns in matrix A
    // A_nnz : Number of nonzero elements in matrix A
    // B : Matrix B
    // B_rownum : Number of rows in matrix B
    // B_colnum : Number of columns in matrix B
    // B_nnz : Number of nonzero elements in matrix B
    // result : return CSR handle for A x B

    MKL_INT job[8];
    job[0] = 0; // the rectangular matrix A is converted to the CSR format;
    job[1] = 0; // zero-based indexing for the rectangular matrix A is used;
    job[2] = 0; // zero-based indexing for the matrix in CSR format is used;
    job[3] = 2; // whole matrix
    //job[4] // maximum number of the non-zero elements allowed if job[0] = 0
    job[5] = 5; // If job[5]>0, arrays acsr, ia, ja are generated for the output storage. If job[5]=0, only array ia is generated for the output storage.
    MKL_INT info = 0; // If info = 0, execution of mkl_zdnscsr was successful.


    MKL_INT status = 1; // return this value to check the execution status  
    //(1 : successfull, 2: error in sparse functions, 3: error in deallocating memory)

    MKL_Complex16 *A_val = (MKL_Complex16 *)mkl_malloc(A_nnz * sizeof(MKL_Complex16), ALIGN);
    MKL_INT *A_col = (MKL_INT *)mkl_malloc(A_nnz * sizeof(MKL_INT), ALIGN);
    MKL_INT *A_row = (MKL_INT *)mkl_malloc( (A_rownum + 1) * sizeof(MKL_INT), ALIGN); // +1 is because we are using 3-aaray variation
    job[4] = A_nnz;
    mkl_zdnscsr(job, &A_rownum, &A_colnum, A, &A_colnum, A_val, A_col, A_row, &info);

    MKL_Complex16 *B_val = (MKL_Complex16 *)mkl_malloc(B_nnz * sizeof(MKL_Complex16), ALIGN);
    MKL_INT *B_col = (MKL_INT *)mkl_malloc(B_nnz * sizeof(MKL_INT), ALIGN);
    MKL_INT *B_row = (MKL_INT *)mkl_malloc((B_rownum + 1) * sizeof(MKL_INT), ALIGN); // +1 is because we are using 3-aaray variation
    job[4] = B_nnz;
    mkl_zdnscsr(job, &B_rownum, &B_colnum, B, &B_colnum, B_val, B_col, B_row, &info);

    sparse_matrix_t   csrA = NULL, csrB = NULL;
    CHECK_SPARSE( mkl_sparse_z_create_csr(&csrA, SPARSE_INDEX_BASE_ZERO, A_rownum, A_colnum, A_row, A_row + 1, A_col, A_val) );
    CHECK_SPARSE( mkl_sparse_z_create_csr(&csrB, SPARSE_INDEX_BASE_ZERO, B_rownum, B_colnum, B_row, B_row + 1, B_col, B_val) );
    CHECK_SPARSE( mkl_sparse_spmm(SPARSE_OPERATION_NON_TRANSPOSE, csrA, csrB, &result) );

memory_free:

    //Release matrix handle and deallocate arrays for which we allocate memory ourselves.
    if (mkl_sparse_destroy(csrA) != SPARSE_STATUS_SUCCESS) status = 3;
    if (mkl_sparse_destroy(csrB) != SPARSE_STATUS_SUCCESS) status = 3;

    //Deallocate arrays for which we allocate memory ourselves.
    mkl_free(A_val); mkl_free(A_col); mkl_free(A_row);
    mkl_free(B_val); mkl_free(B_col); mkl_free(B_row);

    return status;
}

Answer 1

这是工作代码:

/* ***************** Macro ********************* */
#define ALIGN 128

/* To avoid constantly repeating the part of code that checks different functions' status, using the below macros */
#define CHECK_SPARSE(function) do { \
if(function != SPARSE_STATUS_SUCCESS)             \
 {                                                 \
 status = 2;                                       \
 goto memory_free;                                 \
 }                                                 \
} while(0)


/* ****************** Main ******************** */
int main()
{

 << matrices A and B are generated using some data >>

    MKL_INT stat = 0;

    // This part calls the function to multiply matrices as I discussed.

    // A x B --> csrC
    sparse_matrix_t  csrC = NULL;
    stat = dnmm_sp_CSR_handle(CfPrime, Num_of_Buses, Num_of_Branches, CfPrime_nonzero, Yf, Num_of_Branches, Num_of_Buses, Yf_nonzero, &csrC);
    printf("\nstat = %i", stat);

    // Now I convert the csrC to 4-array version of CSR.
    MKL_INT rows, cols;
    sparse_index_base_t indexing = 0;
    MKL_INT *columns_C = NULL, *pointerB_C = NULL, *pointerE_C = NULL;
    MKL_Complex16  *values_C = NULL;
    mkl_sparse_z_export_csr(csrC, &indexing, &rows, &cols, &pointerB_C, &pointerE_C, &columns_C, &values_C);

    // Print the number of rows and columns of converted matrix (which are incorrect sizes)
    printf("\nrows = %i , cols = %i", rows, cols);
}

/* ****************** Function ******************** */
// This function receives two dense matrice, convert them to sparse CSR format, multiply them, and returns the result in CSR handle
int dnmm_sp_CSR_handle(MKL_Complex16 *A, MKL_INT A_rownum, MKL_INT A_colnum, MKL_INT A_nnz, MKL_Complex16 *B, MKL_INT B_rownum, MKL_INT B_colnum, MKL_INT B_nnz, sparse_matrix_t *result) {
// A : Matrix A
// A_rownum : Number of rows in matrix A
// A_colnum : Number of columns in matrix A
// A_nnz : Number of nonzero elements in matrix A
// B : Matrix B
// B_rownum : Number of rows in matrix B
// B_colnum : Number of columns in matrix B
// B_nnz : Number of nonzero elements in matrix B
// result : return CSR handle for A x B

MKL_INT job[8];
job[0] = 0; // the rectangular matrix A is converted to the CSR format;
job[1] = 0; // zero-based indexing for the rectangular matrix A is used;
job[2] = 0; // zero-based indexing for the matrix in CSR format is used;
job[3] = 2; // whole matrix
//job[4] // maximum number of the non-zero elements allowed if job[0] = 0
job[5] = 5; // If job[5]>0, arrays acsr, ia, ja are generated for the output storage. If job[5]=0, only array ia is generated for the output storage.
MKL_INT info = 0; // If info = 0, execution of mkl_zdnscsr was successful.


MKL_INT status = 1; // return this value to check the execution status  
//(1 : successfull, 2: error in sparse functions, 3: error in deallocating memory)

MKL_Complex16 *A_val = (MKL_Complex16 *)mkl_malloc(A_nnz * sizeof(MKL_Complex16), ALIGN);
MKL_INT *A_col = (MKL_INT *)mkl_malloc(A_nnz * sizeof(MKL_INT), ALIGN);
MKL_INT *A_row = (MKL_INT *)mkl_malloc( (A_rownum + 1) * sizeof(MKL_INT), ALIGN); // +1 is because we are using 3-aaray variation
job[4] = A_nnz;
mkl_zdnscsr(job, &A_rownum, &A_colnum, A, &A_colnum, A_val, A_col, A_row, &info);

MKL_Complex16 *B_val = (MKL_Complex16 *)mkl_malloc(B_nnz * sizeof(MKL_Complex16), ALIGN);
MKL_INT *B_col = (MKL_INT *)mkl_malloc(B_nnz * sizeof(MKL_INT), ALIGN);
MKL_INT *B_row = (MKL_INT *)mkl_malloc((B_rownum + 1) * sizeof(MKL_INT), ALIGN); // +1 is because we are using 3-aaray variation
job[4] = B_nnz;
mkl_zdnscsr(job, &B_rownum, &B_colnum, B, &B_colnum, B_val, B_col, B_row, &info);

sparse_matrix_t   csrA = NULL, csrB = NULL;
CHECK_SPARSE( mkl_sparse_z_create_csr(&csrA, SPARSE_INDEX_BASE_ZERO, A_rownum, A_colnum, A_row, A_row + 1, A_col, A_val) );
CHECK_SPARSE( mkl_sparse_z_create_csr(&csrB, SPARSE_INDEX_BASE_ZERO, B_rownum, B_colnum, B_row, B_row + 1, B_col, B_val) );
CHECK_SPARSE( mkl_sparse_spmm(SPARSE_OPERATION_NON_TRANSPOSE, csrA, csrB, result) );

memory_free:

//Release matrix handle and deallocate arrays for which we allocate memory ourselves.
if (mkl_sparse_destroy(csrA) != SPARSE_STATUS_SUCCESS) status = 3;
if (mkl_sparse_destroy(csrB) != SPARSE_STATUS_SUCCESS) status = 3;

//Deallocate arrays for which we allocate memory ourselves.
mkl_free(A_val); mkl_free(A_col); mkl_free(A_row);
mkl_free(B_val); mkl_free(B_col); mkl_free(B_row);

return status;
}