Problem using cusparseScsric0 function

Hi, my name is Ana.

I have been used function in the NVIDIA library: Cusparse and Cublas. I have been implemented in my program the conjugate gradient method using these functions. Thus, it´s worked very well!

However, when I used Cholesky preconditioned and CG iterative methods, in another case, one problem occurred with cusparseScsric0 function. I am using matrix from The University of Florida (Sparse Matrix Collection), website: http://www.cise.ufl.edu/research/sparse/matrices/ in my algorithms with symmetric positive definite linear systems (FIDAP/ex5).

Thus, I have two questions:

1 – Incomplete Cholesky factorization could be: op(A)aproximately R R^T? I read from Matrix Computation´s book (Golub & Van Loan, 2rd Ed., section 10.3.2, pg. 530) about this function.

2 - In the manual “Cuda_cusparse” there is an information about the cusparseScsric0. The function is defined Hermitian/symmetric positive definite sparse matrix (CSR storage format by the three arrays csrValM, csrRowPtrA and csrColIndA). I didn’t understand the supported matrix type: CUSPARSE_MATRIX_TYPE_GENERAL, because only a lower or upper Hermitian/symmetric part of the matrix is actually stored. I used lower or upper or general matrix and I did none result. The results in all the simulations were: CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED.

Please, I need your help!

Answer to 2:
You need to set the Matrix type to CUSPARSE_MATRIX_TYPE_SYMMETRIC ( or CUSPARSE_MATRIX_TYPE_HERMITIAN)
using the routine cusparseSetMatType( cusparseMatDescr_t descrA, cusparseMatrixType_t )

Thank you Philip!!!

I set the parameters as below:

cusparseMatDescr_t descrA = 0;
cusparseStatus = cusparseCreateMatDescr(&descrA); 
cusparseSetMatType(descrA,CUSPARSE_MATRIX_TYPE_SYMMETRIC);
cusparseSetMatIndexBase(descrA,CUSPARSE_INDEX_BASE_ZERO);

In this case, only the main diagonal was calculted!

Do you also defined cusparseSetMatFillMode ? When I define the matrix A as UPPER type, for example, some values were determined incorrectly!

Thanks again!!!

Most of the matrices of the University of Florida are in base 1 ( Fortran legacy )

Could you try :
cusparseSetMatIndexBase(descrA,CUSPARSE_INDEX_BASE_ONE);

For the fillMode, usually they are in LOWER type. But you can easily look into the matrix file to see if it is an UPPER or a LOWER.

Thank you Philip!!!

Another question:

I’m in trouble with the csric0 function!!!

I think there must be a failure on the data structur!!!

CG works very well! But, PCG (cholesky) not converge!!!

I used a small matrix (University Florida - s.p.d.) to check the values. I stored UPPER part of the matrix A for cholesky method. In this case, some values the output function were calculated incorrectly, when I aplly: cusparseScsric0. Some results were: #QNAN0!

Did you use this function correctly?

Thanks, Ana.

Could you please show part of your code?
The matrix is base 1 and lower triangle, but you can use csr2coo to have a upper triangle one.

Hi,

Same things for me. I tried with a symmetric matrix where I store only the upper side. The base index is 0 (I made it myself). Some values are indefinite after using csric0.
There is something strange: I store the upper side of the matrix but if I say in the descriptor that the matrix is lower I don’t have anymore indefinite value. Both cases, the PCG does not converge.

EDIT: Actualy, if I put in the descriptor that the matrix is lower instead of upper, it returns me exactly the same matrix.
I guess, there is a problem when I use csric0 for the upper side matrix.

Seltymar, could you show your code? Otherwise it is not easy to find the root cause.

This is the code used to compute the incomplete cholesky.
The matrix is composed of dev_valR, dev_row_ptrR and dev_colIndR which correspond to the CSR matrix. I tested the matrix on a PCG using a Jacobi preconditioner. It is slow but works.
Also, it happens often that using cusparseScsrsv_analysis and cusparseScsric0 makes the cusparseScsrsv_analysis fail with the code error 6: execution failed and the cusparseScsric0 has the error 7: internal error.

int N = width*height;

/* Get handle to the CUBLAS context */
    cublasHandle_t cublasHandle = 0;
    cublasStatus_t cublasStatus;
    cublasStatus = cublasCreate(&cublasHandle);
    if(cublasStatus != CUBLAS_STATUS_SUCCESS)
        fprintf(stderr, "cublasCreate returned error code %d !
", cublasStatus);

    /* 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 !
", cusparseStatus);

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

    cusparseSetMatFillMode(descrR,CUSPARSE_FILL_MODE_UPPER);
    cusparseSetMatType(descrR,CUSPARSE_MATRIX_TYPE_SYMMETRIC);
    cusparseSetMatIndexBase(descrR,CUSPARSE_INDEX_BASE_ZERO);

    // create the analysis info object for the R matrix 
    cusparseSolveAnalysisInfo_t infoR = 0;
    cusparseStatus = cusparseCreateSolveAnalysisInfo(&infoR);
    if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
        fprintf(stderr, "cusparseCreateSolveAnalysisInfo returned error code %d !
", cusparseStatus);
	
	
    cusparseStatus = cusparseScsrsv_analysis(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, nnz_Symmetric, descrR, dev_valR, dev_row_ptrR, dev_colIndR, infoR);
    if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
        fprintf(stderr, "cusparseScsrsv_analysis returned error code %d !
", cusparseStatus);
	
	
    // generate the Incomplete Cholesky factor H for the matrix R using cusparseScsric0 
    cusparseStatus = cusparseScsric0(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, descrR, dev_valR, dev_row_ptrR, dev_colIndR, infoR);
    if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
        fprintf(stderr, "cusparseScsric0 returned error code %d !
", cusparseStatus);

The code looks fine.
If csrsv_analysis failed, either csr format is wrong or there exists some bug.
If (dev_rowPtrR, dev_colIndR, dev_valR) is indeed upper triangular matrix of dimension NxN, nnz_Symmetric nonzero entries, we need to take a look on the matrix.
Could you please file a bug?
(If you are not a registered developer, please sign up, then file a bug.)

As the cusparse function was correct, I focused and checked again my matrix and I fixed a bug so now the incomplete cholesky looks correct (no more indefinite values).
Also, in the documentation it is written that the dev_rowPtrR has m+1 which contains the beginning of all row and the end of the last one+1. My pointer contained only the beginning of each row, I changed it to have the last element.
But I still have some problem with my residual. I will check it again.
This is the rest of the conjugate gradient (I didn’t follow exactly what Maxim Naumov did):

// create the info and analyse the lower and upper triangular factors
cusparseSolveAnalysisInfo_t inforRt = 0;
cusparseSolveAnalysisInfo_t inforR = 0;
cusparseCreateSolveAnalysisInfo(&inforRt ) ;
cusparseCreateSolveAnalysisInfo(&inforR ) ;
cusparseStatus=cusparseScsrsv_analysis(cusparseHandle , CUSPARSE_OPERATION_TRANSPOSE ,N , nnz_Symmetric, descrR , dev_valR , dev_row_ptrR , dev_colIndR , inforRt);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
    fprintf(stderr, "cusparseScsrsv_analysis inforRt returned error code %d !", cusparseStatus);

cusparseStatus=cusparseScsrsv_analysis(cusparseHandle , CUSPARSE_OPERATION_NON_TRANSPOSE ,N , nnz_Symmetric, descrR , dev_valR , dev_row_ptrR , dev_colIndR , inforR);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
    fprintf(stderr, "cusparseScsrsv_analysis inforR returned error code %d !", cusparseStatus);
	
//  Ax = A*x
cusparseStatus = cusparseScsrmv(cusparseHandle,CUSPARSE_OPERATION_NON_TRANSPOSE, N, N, nnz, &alpha, descr, dev_val, dev_row_ptr, dev_colInd, dev_alphaLaplacian, &beta, dev_Ax);

// r= r - Ax
cublasStatus = cublasSaxpy(cublasHandle, N, &alpham1, dev_Ax, 1, d_r, 1);

// z_{0} = C^{-1}*r_{0}
cusparseScsrsv_solve(cusparseHandle, CUSPARSE_OPERATION_TRANSPOSE , N, &alpha, descrR, dev_valR, dev_row_ptrR, dev_colIndR, inforRt, d_r, d_y);
cusparseScsrsv_solve(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE , N,   &alpha, descrR, dev_valR, dev_row_ptrR, dev_colIndR, inforR, d_y, dev_z);

cublasStatus = cublasSdot(cublasHandle, N, d_r, 1, dev_z, 1, &r1);

k = 1;
	
printf("residual = %e", sqrt(r1));
	
while (r1 > tol*tol && k <= 20 )//max_iter)
{	
    if (k > 1) {
        //b_{i} = (r_{i+1}, z_{i+1}/(r_{i} , z_{i})
        b = r1 / r0;

        // p_{i+1}= z_{i+1} + b_{i}* p_{i} 
        cublasStatus = cublasSscal(cublasHandle, N, &b, dev_p, 1);

        cublasStatus = cublasSaxpy(cublasHandle, N, &alpha, dev_z, 1, dev_p, 1);
    }  else {
        // p = z
        cublasStatus = cublasScopy(cublasHandle, N, dev_z, 1, dev_p, 1);
    }
		
    // Ax_{i} = A*p^{T}_{i}
    cusparseScsrmv(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, N, nnz, &alpha, descr, dev_val, dev_row_ptr, dev_colInd, dev_p, &beta, dev_Ax);

     // Ax_{i} = p_{i}*Ax_{i}
     cublasStatus = cublasSdot(cublasHandle, N, dev_p, 1, dev_Ax, 1, &dot);

		
      //a = (r_{i} , z_{i})/(Ap_{i} , p_{i})
      a = r1 / dot;
		
      // x_{i+1}= x_{i} + a*p_{i}
      cublasStatus = cublasSaxpy(cublasHandle, N, &a, dev_p, 1, dev_alphaLaplacian, 1);

na = -a;

	// r_{i+1} = r_{i} - a *Ax_{i}
        cublasStatus = cublasSaxpy(cublasHandle, N, &na, dev_Ax, 1, d_r, 1);

        r0 = r1;

    // z_{i+1} = C^{-1}*r_{i+1}		
    cusparseScsrsv_solve(cusparseHandle, CUSPARSE_OPERATION_TRANSPOSE , N, &alpha, descrR, dev_valR, dev_row_ptrR,    dev_colIndR,inforRt, d_r, d_y);
     cusparseScsrsv_solve(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE , N,   &alpha, descrR, dev_valR, dev_row_ptrR, dev_colIndR, inforR, d_y, dev_z);
	
      // r1_{i+1} = r_{i+1} * z_{i+1}
      cublasStatus = cublasSdot(cublasHandle, N, d_r, 1, dev_z, 1, &r1);
        printf("iteration = %3d, residual = %e", k, sqrt(r1));	
	k++;		
}

The matrix (dev_val, dev_row_ptr, dev_colInd) is the same matrix than dev_valR but containing all data. I will try the method of maxim Naumov, maybe it will resolve my problem.

EDIT:
My residual r1 is < 0, so I guess my matrix R is still not correct.

Seltymar, I think that has one problem with: cusparseScsrsv_solve and cusparseScsric0. My program not converged!!!

I resolved CG and its ok!

When aplly PCG not converged.
I calculated incomplete cholesky separately (not use cusparseScsric0) and aplly cusparseScsrsv_solve (preconditioner) function and not converged (same matrix resolved in CG). In this case, PCG, I added just cusparseScsrsv_solve function. Thus, there is one problem!!!

When aplly cusparseScsric0 same things for me!!!

Ana.

Hi Ana,
Do you still have indefinite value like #QNANO after using cusparseScsric0 ? I have correct values now (My matrix was incorrect but converged with PCG using a Jacobi preconditioner so it’s difficult to know what’s going wrong).

I will check my result after using cusparseScrsv_solve and if I find nothing I will file a bug.

EDIT:

cusparseScsrsv_solve(cusparseHandle, CUSPARSE_OPERATION_TRANSPOSE , N, &alpha, descrR, dev_valR, dev_row_ptrR, dev_colIndR, inforRt, d_r, d_y);

The result of this function is not correct. I have wrong values in the output d_y.

I applied to the CUDA developper program. I’m waiting to be accepted.

Hi Seltymar,

When aplly cusparseScsric0, I found some values with: #QNANO.
My matrix is small (27 x 27) and I have wrong values in the output this function.

Set:

cusparseSetMatIndexBase(descrA,CUSPARSE_INDEX_BASE_ZERO);
cusparseSetMatType(descrA,CUSPARSE_MATRIX_TYPE_SYMMETRIC);
cusparseSetMatFillMode(descrA,CUSPARSE_FILL_MODE_UPPER);

Did you program converged, when used cusparseScsric0 whit CG?

, Ana.

I have correct values when using cusparseScsrcric0 but it does not converge.
For me, the function cusparseScrsv_solve does nothing. If I initialize my output vector to 5 for example, after using cusparseScrsv_solve, it gives me 5.

I have same parameter for the matrix.

EDIT:
I have filed a bug about cusparseScrsv_solve.

EDIT 2:
I’m not sure if it is interesting to use PCG with a incomplete Cholesky. The time for the analyze and to create the Cholesky matrix takes around the same time than a PCG using a Jacobi preconditioner which is fast to invert.
To execute this

// z_{0} = C^{-1}*r_{0}

, I created a kernel function with CUDA and it’s quite fast. But of course, it depends how many iterations you CG does. I have around 150 iterations which is not so big.

In my case, cusparseScsrcric0 function is with some wrong values.
How did you used this function (parameter)? Did you change something in this function?

I configured:

cusparseSetMatIndexBase(descrA,CUSPARSE_INDEX_BASE_ZERO);
cusparseSetMatType(descrA,CUSPARSE_MATRIX_TYPE_SYMMETRIC);
cusparseSetMatFillMode(descrA,CUSPARSE_FILL_MODE_UPPER);

I think that using cusparseSetMatFillMode like UPPER appears a problem. How did you used?

Thanks, Ana.

i have set parameters like this:

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

cusparseSetMatFillMode(descrR,CUSPARSE_FILL_MODE_UPPER);
cusparseSetMatType(descrR,CUSPARSE_MATRIX_TYPE_SYMMETRIC);
cusparseSetMatIndexBase(descrR,CUSPARSE_INDEX_BASE_ZERO);
 
// create the analysis info object for the R matrix 
cusparseSolveAnalysisInfo_t infoR = 0;
cusparseStatus = cusparseCreateSolveAnalysisInfo(&infoR);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
     fprintf(stderr, "cusparseCreateSolveAnalysisInfo returned error code %d !
", cusparseStatus);
	
cusparseStatus = cusparseScsrsv_analysis(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, nnz_Symmetric, descrR, dev_valR, dev_row_ptrR, dev_colIndR, infoR);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
    fprintf(stderr, "cusparseScsrsv_analysis returned error code %d !
", cusparseStatus);
	
// generate the Incomplete Cholesky factor H for the matrix R using cusparseScsric0 
cusparseStatus = cusparseScsric0(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, N, descrR, dev_valR, dev_row_ptrR, dev_colIndR, infoR);
if(cusparseStatus!=CUSPARSE_STATUS_SUCCESS)
    fprintf(stderr, "cusparseScsric0 returned error code %d !
", cusparseStatus);

the matrix R (csr format: dev_valR, dev_row_ptrR, dev_colIndR) is size N with nnz_Symmetric non zero elements. The matrix R is upper side of the symmetric matrix.

Its like my parameters, except dev_row_ptrR array where I put (N+1) values. Did your array is N?

How did you set nnz_Symmetric? Its just UPPER ou full symmetric matrix.

Your matrix is small? I´m testing 27 x 27 matrix with few nonzeros values.

cusparseSolveAnalysisInfo_t infoA = 0;
cusparseStatus = cusparseCreateSolveAnalysisInfo(&infoA);

cusparseMatDescr_t descrA = 0;
cusparseStatus = cusparseCreateMatDescr(&descrA);
cusparseSetMatIndexBase(descrA,CUSPARSE_INDEX_BASE_ZERO);
cusparseSetMatType(descrA,CUSPARSE_MATRIX_TYPE_SYMMETRIC);
cusparseSetMatFillMode(descrA,CUSPARSE_FILL_MODE_UPPER);

cusparseStatus = cusparseScsrsv_analysis(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, M, nnz, descrA, d_ValAICP, d_csrRowPtrA, d_cooColIndA, infoA);
if ( checkCusparseStatus (cusparseStatus, (char*)“cusparseScsrv_analysis” ) ) return EXIT_FAILURE;

cusparseStatus = cusparseScsric0(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, M, descrA, d_ValAICP, d_csrRowPtrA, d_cooColIndA, infoA);
if ( checkCusparseStatus (cusparseStatus, (char*)“cusparseScsric0” ) ) return EXIT_FAILURE;

Thanks, Ana.

I have (N+1) values for the dev_roz_ptrR.

nnz_Symmetric is the number of non-zero values only for the upper side.

I have test it with several matrices but they are all big, several millions of nonzeros values.

I have the same code than you, if nnz is the number of nonzeros values only for the upper side.

Maybe if you can give me your input data I can try myself.

Seltymar,

The matrix that I´m using: http://www.cise.ufl.edu/research/sparse/matrices/FIDAP/ex5.html

I inverted rows and columns to obtain UPPER side.

Thanks, I await your response!!!

Ana.