Search code examples
cudacublas

cusparse csrsvanalysis sometimes works, sometimes fails

I'm trying to use the preconditoned conjugate gradient to resolve Ax=b. So I took example on the sample given with the cuda-sdk. Sometimes, when I call the function cusparseScsrsv_analysis, it returns the error 6 which is "execution failed". Sometimes, it works.

The matrix A is symmetric positive definite.

Also, the conjugate gradient works correctly on the same data.

Here is my code: 

/* Get handle to the CUSPARSE context */
cusparseHandle_t cusparseHandle = 0;
cusparseStatus_t cusparseStatus;
cusparseStatus = cusparseCreate(&cusparseHandle);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
    fprintf(stderr, "cusparseCreate returned error code %d !\n", cusparseStatus);

cusparseMatDescr_t descr = 0;
cusparseStatus = cusparseCreateMatDescr(&descr);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
    fprintf(stderr, "cusparseCreateMatDescr returned error code %d !\n", cusparseStatus);

// create the analysis info object for the A matrix 
cusparseSolveAnalysisInfo_t infoA = 0;
cusparseStatus = cusparseCreateSolveAnalysisInfo(&infoA);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
    fprintf(stderr, "cusparseCreateSolveAnalysisInfo returned error code %d !\n", cusparseStatus);

// Perform the analysis for the Non-Transpose case 
cusparseStatus = cusparseScsrsv_analysis(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, nnz, descr, dev_val, dev_row_ptr, dev_colInd, infoA);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
    fprintf(stderr, "cusparseScsrsv_analysis 1 returned error code %d !\n", cusparseStatus);

N is the number of column and row, nnz is the number of non-zero elements. My matrix is in csr format.

EDIT: I don't see any special requirement. I don't think this this is du to memory, i have more than 2GB and I'm not using a big matrix (48MB).

I tried the preconditionned conjugate gradient with a jacobi precondionner and it works correctly too but if I try to analyze with cusparse, it fails half the time.

What I want is to use the algorithm of Maxim Noumov (http://developer.download.nvidia.com/compute/DevZone/docs/html/CUDALibraries/doc/Preconditioned_Iterative_Methods_White_Paper.pdf) using cusparse and cublas.

EDIT2:

I need some explanation about curspace. If I put in the descriptor this line: cusparseSetMatType(descr,CUSPARSE_MATRIX_TYPE_SYMMETRIC); the analyse works but the weird thing is that I store the whole matrix not only the upper or lower part. If I put cusparseSetMatType(descr,CUSPARSE_MATRIX_TYPE_GENERAL);, it doesn't work. Moreover I don't understand why I have to store in the dev_row_ptr m+1 elements where m is the number of row. What do I put in the last element ?

Other question: the function cusparseScsric0 takes as input/output the matrix value （csrValM in the documentation) which is the whole matrix as input and the incomplete-CHolesky upper or lower triangular only as output. How does it work ?

Solution

  • The cusparse documentation about cusparseScsric0 is wrong, it takes CUSPARSE_MATRIX_TYPE_SYMMETRIC as input. This function made the cusparseScsrsv_analysis crashed.

    Here is a correct code: 

    cusparseMatDescr_t descrR = 0;
cusparseStatus = cusparseCreateMatDescr(&descrR);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
fprintf(stderr, "cusparseCreateMatDescr returned error code %d !\n", cusparseStatus);

cusparseSetMatFillMode(descrR,CUSPARSE_FILL_MODE_UPPER); // It can also be lower side
cusparseSetMatType(descrR,CUSPARSE_MATRIX_TYPE_SYMMETRIC);
cusparseSetMatIndexBase(descrR,CUSPARSE_INDEX_BASE_ZERO);

cusparseSolveAnalysisInfo_t infoR = 0;
cusparseStatus = cusparseCreateSolveAnalysisInfo(&infoR);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
    fprintf(stderr, "cusparseCreateSolveAnalysisInfo returned error code %d !\n", cusparseStatus);


cusparseStatus = cusparseScsrsv_analysis(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, 27, 153, descrR, dev_valR, dev_row_ptrR, dev_colIndR, infoR);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
    fprintf(stderr, "cusparseScsrsv_analysis returned error code %d !\n", cusparseStatus);

// generate the Incomplete Cholesky factor H for the matrix R using cusparseScsric0 
cusparseStatus = cusparseScsric0(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, 27, descrR, dev_valR, dev_row_ptrR, dev_colIndR, infoR);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
    fprintf(stderr, "cusparseScsric0 returned error code %d !\n", cusparseStatus);

    Also dev_row_ptrR is number of row+1.