除非我增加稀疏矩阵的稀疏度，否则 Cusparse 非法内存访问

Question

我正在尝试制作一个现有的软件，该软件使用手动调整的特殊 CSC 矩阵的稀疏乘法，每列恰好有 k 个非零元素。我决定使用 cusparse 来完成这项工作，但不幸的是我遇到矩阵乘法在某些情况下需要超过 7 秒，这比代码的 CPU 版本慢得多。 （最大稀疏矩阵为 19871x1000，最大密集矩阵为 1000*150，nnz = 101000）。

当尝试在一个自包含示例中重现该问题时，我总是在 nnz != sparse_cols 时遇到“非法内存访问错误”。

经过一番调查后发现，如果我将矩阵的大小增加 10 倍，问题就会消失。如果我使矩阵足够小，我就不会遇到崩溃。然而，对于大型矩阵，稀疏矩阵必须不跨越某种程度的密集度，否则乘法会产生一堆非法内存访问。 这是出现问题的代码：

#include <cuda.h>
#include <cusparse.h>
#include <iostream>
#include <stdlib.h> 

#define CALL_CUDA( err ) \
{ if (err != cudaSuccess) \
    {std::cout<<"cuda Error "<< cudaGetErrorString(err)<<" in "<<__FILE__<<" at line "<<__LINE__<<"\n"; exit(EXIT_FAILURE); }\
}


int main(){

    //cusparse status and handle
    cusparseStatus_t status; 
    cusparseHandle_t handle = 0;
    status = cusparseCreate(&handle);

    if (status != CUSPARSE_STATUS_SUCCESS){
        std::cout << "Error creating handle: " << status << std::endl;
    }

    //Set matrix description
    cusparseMatDescr_t descr; //Describe the matrices
    cusparseCreateMatDescr(&descr);
    cusparseSetMatType(descr, CUSPARSE_MATRIX_TYPE_GENERAL);
    cusparseSetMatIndexBase(descr, CUSPARSE_INDEX_BASE_ZERO);

    //Sparse matrix properties
    int sparse_rows = 19871;
    int sparse_cols = 1000;
    int nnz_new = 101000;
    //int nnz_new = 1000; //Works with that value

    //Dense matrix properties
    int bmat_rows = 1000;
    int bmat_cols = 150;


    //Generate a special type of sparse matrix that has exactly k nonzero elements in each column in CSC format
    float * amat_vals;
    CALL_CUDA(cudaMallocHost((void **)&amat_vals, nnz_new*sizeof(float)));
    int * amat_idx;
    CALL_CUDA(cudaMallocHost((void **)&amat_idx, nnz_new*sizeof(int)));
    int * crccolptr;
    CALL_CUDA(cudaMallocHost((void **)&crccolptr, (sparse_cols+1)*sizeof(int)));

    //Fill in values with random values
    for (int i = 0; i < nnz_new; i++){
        amat_vals[i] = (float)rand()/(float)RAND_MAX;
    }

    //Generate indexes for those rows
    for (int i = 0; i < nnz_new; i++){
        amat_idx[i] = rand() % (sparse_rows - 1);
    }

    //generate crccolptr
    int k = (int)(nnz_new/sparse_cols); //Number of elements per row
    for (int i = 0; i < sparse_cols; i++){
        crccolptr[i] = k*i;
    }
        crccolptr[sparse_cols] = nnz_new;


    //Generate bmat_array with random floats
    float * bmat_array;
    CALL_CUDA(cudaMallocHost((void **)&bmat_array, bmat_rows*bmat_cols*sizeof(float))); 
    for (int i = 0; i < bmat_rows*bmat_cols; i++){
        bmat_array[i] = (float)rand()/(float)RAND_MAX;
    }

    //generate array for output
    float * outmatrix_test;
    CALL_CUDA(cudaMallocHost((void **)&outmatrix_test, sparse_rows*bmat_cols*sizeof(float)));

    //Allocate and copy device memory for sparse matrix
    float * cudavals;
    int * colIdx;
    int * colPtr;

    CALL_CUDA(cudaMalloc((void **)&colPtr, (sparse_cols + 1)*sizeof(int)));
    CALL_CUDA(cudaMemcpy(colPtr, crccolptr, (sparse_cols + 1)*sizeof(int), cudaMemcpyHostToDevice));

    CALL_CUDA(cudaMalloc((void **)&cudavals, nnz_new*sizeof(float)));
    CALL_CUDA(cudaMalloc((void **)&colIdx, nnz_new*sizeof(int)));

    CALL_CUDA(cudaMemcpy(cudavals, amat_vals, nnz_new*sizeof(float), cudaMemcpyHostToDevice));
    CALL_CUDA(cudaMemcpy(colIdx, amat_idx, nnz_new*sizeof(int), cudaMemcpyHostToDevice));

    //Allocate and copy device memory for dense matrix
    float * B_gpumatrix;
    CALL_CUDA(cudaMalloc((void **)&B_gpumatrix, bmat_rows*bmat_cols*sizeof(float)));
    CALL_CUDA(cudaMemcpy(B_gpumatrix, bmat_array, bmat_rows*bmat_cols*sizeof(float), cudaMemcpyHostToDevice));

    //Allocate output matrix
    float * outmatrix_gpu; 
    CALL_CUDA(cudaMalloc((void **)&outmatrix_gpu, (sparse_rows*bmat_cols)*sizeof(float)));

    //sparse_cols is passed as sparse_rows, because we're multiplying a CSC matrix instead of a CSR so we need
    // to transpose it and invert the rows and columns.
    const float alpha = 1.0;
    const float beta = 0.0;

/*
    float * outmatrix_gpu2; 
    CALL_CUDA(cudaMalloc((void **)&outmatrix_gpu2, (sparse_rows*sparse_cols)*sizeof(float)));
    cusparseStatus_t mat_mul = cusparseScsc2dense(handle, sparse_rows, sparse_cols, descr, cudavals, colIdx, colPtr, outmatrix_gpu2, sparse_rows);
    float * outmatrix_test2;
    CALL_CUDA(cudaMallocHost((void **)&outmatrix_test2, sparse_rows*sparse_cols*sizeof(float)));
    CALL_CUDA(cudaMemcpy(outmatrix_test2, outmatrix_gpu2, (sparse_rows*sparse_cols)*sizeof(float), cudaMemcpyDeviceToHost));
*/
    cusparseStatus_t mat_mul = cusparseScsrmm(handle,  //Cusparse handle
     CUSPARSE_OPERATION_TRANSPOSE, //Transposing the matrix
     sparse_cols, //Number of sparse rows. Since we're using CSC matrix it's the columns.
     bmat_cols,   //Number of columns of the dense matrix
     sparse_rows, //Number of sparse cols, Since we're using CSC matrix it's the rows
     nnz_new,     //Non zero elements
     &alpha,      //Pointer to alpha (1.0)
     descr,       //Description of the matrix
     cudavals,    //The values vector
     colPtr,      //The column pointer
     colIdx,      //The indexes of the sparse matrix
     B_gpumatrix, //Dense matrix array
     bmat_rows,   //ldb - the rows of the dense matrix
     &beta,       //Pointer to beta. It's 0
     outmatrix_gpu, //Pointer to the output matrix
     sparse_rows); //ldc - leading dimensions of the output matrix.

    if (mat_mul != CUSPARSE_STATUS_SUCCESS){
        std::cout << "MULTIPLICATION ERROR: " << mat_mul << std::endl;
    }
    cudaThreadSynchronize(); //Syncs before copy back to memory should not be necessary
    cudaDeviceSynchronize();

    //Copy matrix back to host
    CALL_CUDA(cudaMemcpy(outmatrix_test, outmatrix_gpu, (sparse_rows*bmat_cols)*sizeof(float), cudaMemcpyDeviceToHost));

    CALL_CUDA(cudaFree(outmatrix_gpu));
    CALL_CUDA(cudaFree(cudavals));
    CALL_CUDA(cudaFree(colPtr));
    CALL_CUDA(cudaFree(colIdx));
    CALL_CUDA(cudaFree(B_gpumatrix));
    CALL_CUDA(cudaFreeHost(crccolptr));
    CALL_CUDA(cudaFreeHost(amat_vals));
    CALL_CUDA(cudaFreeHost(amat_idx));
    CALL_CUDA(cudaFreeHost(bmat_array));
    CALL_CUDA(cudaFreeHost(outmatrix_test));

    return 1;
}

我相信我正在生成一个有效的稀疏矩阵，因为我可以使用适当的 cusparse 函数将其转换为密集矩阵，而不会触发任何无效的内存访问。

在 cuda-memcheck 下运行上述代码时，您可以看到来自 cusparseScsrmm 的许多非法访问。在没有 cuda-memcheck 的情况下运行，您会在矩阵乘法后的第一个 cuda 操作中看到错误。

任何想法我做错了什么？我希望如果我能解决这个问题，我将能够诊断（或至少隔离）一个自包含的示例，该示例展示了令人痛苦的缓慢的矩阵乘法。

编辑：

使用较小的矩阵我没有遇到问题。具有 50*200 的稀疏矩阵对 NNZ 工作正常，直到大约 1000，但在 NNZ = 5000 时需要永远（我在半分钟后将其杀死）。将矩阵大小增加到 200*500 可以在 NNZ = 5000 的情况下立即执行......奇怪。

EDIT2：

如果我将矩阵的大小增加 10 倍，则 nnz 的原始数量有效。

Answer 1

这是不明智的：

//Generate indexes for those rows
for (int i = 0; i < nnz_new; i++){
    amat_idx[i] = rand() % (sparse_rows - 1);
}

CSR matrix format expects 要按从左到右、从上到下的顺序存储的值向量。因此，每行中的列索引必须按升序。您正在以随机顺序生成列索引，事实上，您很可能会在同一行中生成两个元素具有相同的列索引。那简直是坏掉了。

您的变量命名也让我有些困惑。 CSR 是压缩的稀疏行格式，它期望：

1. 矩阵值向量（长度=nnz）
2. 列索引的向量指定每个值属于哪一列（=nnz 长度）
3. 行指针的向量指定每行的开始（=numrows +1 长度）

由于您使用的是Scsrmm 函数，CSR 格式是必需的。

像 crccolptr 这样的变量名称在 CSR 格式中对我来说没有意义。

作为一个简单的证明点，将上面摘录的代码替换为以下内容：

//Generate indexes for those rows
int my_idx = 0;
int j;
for (int i = 0; i < sparse_rows; i++){
    //amat_idx[i] = rand() % (sparse_rows - 1);
    for (j = 0; j < (nnz_new/sparse_rows); j++)
      amat_idx[my_idx++] = j;

}
while (my_idx < nnz_new) amat_idx[my_idx++] = j++;

而且我相信错误会消失，因为实际矩阵现在符合 CSR 格式的预期。