cusparse csrsvanalysis fonctionne parfois, échoue parfois

Question

J'essaie d'utiliser le gradient conjugué préconditoné pour résoudre Ax = b. J'ai donc pris exemple sur l'exemple donné avec le cuda-sdk. Parfois, quand j'appelle la fonction cusparseScsrsv_analysis, il renvoie l'erreur 6 qui est "exécution échouée". Parfois, cela fonctionne.

La matrice A est définie positive symétrique.

De plus, le gradient conjugué fonctionne correctement sur les mêmes données.

Voici mon 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 est le nombre de colonnes et de rangées, nnz est le nombre d'éléments non nuls. Ma matrice est au format CSR.

EDIT: Je ne vois aucune exigence particulière. Je ne pense pas que cela soit du à la mémoire, j'ai plus de 2Go et je n'utilise pas une grosse matrice (48MB).

J'ai essayé le gradient conjugué préconditionné avec un préconditionneur jacobi et cela fonctionne correctement aussi, mais si j'essaie d'analyser avec cusparse, il échoue la moitié du temps. Ce que je veux, c'est utiliser l'algorithme de Maxim Noumov (http://developer.download.nvidia.com/compute/DevZone/docs/html/CUDALibraries/doc/Preconditioned_Iterative_Methods_White_Paper.pdf) en utilisant cusparse et cublas.

EDIT2:

J'ai besoin d'explications sur curspace. Si je mets dans le descripteur cette ligne: cusparseSetMatType(descr,CUSPARSE_MATRIX_TYPE_SYMMETRIC); l'analyse fonctionne mais la chose étrange est que je stocke la matrice entière non seulement la partie supérieure ou inférieure. Si je mets cusparseSetMatType(descr,CUSPARSE_MATRIX_TYPE_GENERAL);, cela ne fonctionne pas. En outre, je ne comprends pas pourquoi je dois stocker dans le dev_row_ptr m + 1 éléments où m est le numéro de ligne. Qu'est-ce que je mets dans le dernier élément?

Autre question: la fonction cusparseScsric0 prend en entrée/sortie la valeur de la matrice (csrValM dans la documentation) qui est la matrice entière en entrée et la triangulaire supérieure ou inférieure incomplète-Cholesky uniquement en tant que sortie. Comment ça marche ?

Answer 1

La documentation de cusparse sur cusparseScsric0 est erronée, il prend CUSPARSE_MATRIX_TYPE_SYMMETRIC en entrée. Cette fonction a causé la panne du cusparseScsrsv_analysis.

Voici un code correct:

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); 

également dev_row_ptrR est le numéro de la ligne + 1.