【发布时间】：2018-02-16 01:58:41
【问题描述】：
我正在尝试开发一个将矩阵 A 和 B 相乘的函数，它们的格式一般但本质上是稀疏的。这些矩阵包含复数。我的问题是，当我不使用该函数并在 main() 中编写所有内容时，乘法对于任何大小的数组都非常有效。但是当我使用自己的函数时，结果是损坏的，并且大多数时候我会收到随机错误消息。
这是这个函数的作用：
- 将 A 和 B 转换为 CSR 格式 (mkl_zdnscsr)。
- 使用步骤 1 中的数据为 A 和 B 创建 CSR 句柄 (mkl_sparse_z_create_csr)。
- 在步骤 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;
}
【问题讨论】：
-
我现在可以确认问题了。这是从函数返回“结果”的方式。如果我使用 mkl_sparse_z_export_csr 导出函数内部“结果”的 CSR 向量（在 mkl_sparse_spmm 之后），一切正常。但是，如果我在 main() 内部和调用函数后立即做同样的想法，我会得到损坏的数据。