5.2 CUSPARSE

CUSPARSE — библиотека для CUDA, содержащая набор базовых подпрограмм линейной алгебры, используемых для обработки разреженных матриц, которая предназначена для вызова из C или C++. Эти подпрограммы могут быть разделены на четыре категории:

1.
Подпрограммы 1-го уровня включают в себя операции между разреженными векторами и векторами в плотной форме.
2.
Подпрограммы 2-го уровня включают в себя операции между разреженными матрицами и векторами в плотной форме.
3.
Подпрограммы 3-го уровня включают в себя операции между разреженными матрицами и множеством векторов-столбцов в плотной форме.
4.
Преобразовательные подпрограммы, которые позволяют преобразовывать матрицы в различные форматы.

Библиотека написана с использованием модели параллельного программирования CUDA и использует преимущество вычислительных ресурсов графических процессоров NVIDIA (GPU). Интерфейс программирования приложений CUSPARSE предполагает, что входные и выходные данные находятся в памяти GPU (устройство), а не в памяти процессора (хост), до тех пор пока память процессора точно не указана в строке HostPtr, являясь частью имени параметра функции. Выделение памяти и обмен данными между памятью GPU и памятью процессора с помощью стандартных подпрограмм CUDA, таких как, cudaMalloc(), cudaFree(), cudaMemcpy(), и cudaMemcpyAsync() осуществляет сам пользователь. Библиотека в настоящее время работает только на системах с одной видеокартой, и не распараллеливается между несколькими GPU.

Пример, который отображает основные возможности библиотеки CUSPARSE:

#include <stdio.h> 
#include <stdlib.h> 
#include <cuda_runtime.h> 
#include "CUSPARSE.h" 
#define CLEANUP(s) \ 
do { \ 
printf ("%s\n", s); \ 
if (yHostPtr) free(yHostPtr); \ 
if (zHostPtr) free(zHostPtr); \ 
if (xIndHostPtr) free(xIndHostPtr); \ 
if (xValHostPtr) free(xValHostPtr); \ 
if (cooRowIndexHostPtr) free(cooRowIndexHostPtr);\ 
if (cooColIndexHostPtr) free(cooColIndexHostPtr);\ 
if (cooValHostPtr) free(cooValHostPtr); \ 
if (y) cudaFree(y); \ 
if (z) cudaFree(z); \ 
if (xInd) cudaFree(xInd); \ 
if (xVal) cudaFree(xVal); \ 
if (csrRowPtr) cudaFree(csrRowPtr); \ 
if (cooRowIndex) cudaFree(cooRowIndex); \ 
if (cooColIndex) cudaFree(cooColIndex); \ 
if (cooVal) cudaFree(cooVal); \ 
if (handle) CUSPARSEDestroy(handle); \ 
fflush (stdout); \ 
} while (0) 
int main(){ 
cudaError_t cudaStat1,cudaStat2,cudaStat3,cudaStat4,cudaStat5,cudaStat6; 
CUSPARSEStatus_t status; 
CUSPARSEHandle_t handle=0; 
CUSPARSEMatDescr_t descra=0; 
int * cooRowIndexHostPtr=0; 
int * cooColIndexHostPtr=0; 
double * cooValHostPtr=0; 
int * cooRowIndex=0; 
int * cooColIndex=0; 
double * cooVal=0; 
int * xIndHostPtr=0; 
double * xValHostPtr=0; 
double * yHostPtr=0; 
int * xInd=0; 
double * xVal=0; 
double * y=0; 
int * csrRowPtr=0; 
double * zHostPtr=0; 
double * z=0; 
int n, nnz, nnz_vector, i, j; 
printf("testing example\n");

Cозданим разреженную матрицу в формате COO:

/* |1.0     2.0 3.0| 
   |    4.0        | 
   |5.0     6.0 7.0| 
   |    8.0     9.0| */ 
n=4; nnz=9; 
cooRowIndexHostPtr = (int *) malloc(nnz*sizeof(cooRowIndexHostPtr[0])); 
cooColIndexHostPtr = (int *) malloc(nnz*sizeof(cooColIndexHostPtr[0])); 
cooValHostPtr = (double *)malloc(nnz*sizeof(cooValHostPtr[0])); 
if ((!cooRowIndexHostPtr) || (!cooColIndexHostPtr) || (!cooValHostPtr)){ 
CLEANUP("Host malloc failed (matrix)"); 
return EXIT_FAILURE; 
} 
cooRowIndexHostPtr[0]=0; cooColIndexHostPtr[0]=0; cooValHostPtr[0]=1.0; 
cooRowIndexHostPtr[1]=0; cooColIndexHostPtr[1]=2; cooValHostPtr[1]=2.0; 
cooRowIndexHostPtr[2]=0; cooColIndexHostPtr[2]=3; cooValHostPtr[2]=3.0; 
cooRowIndexHostPtr[3]=1; cooColIndexHostPtr[3]=1; cooValHostPtr[3]=4.0; 
cooRowIndexHostPtr[4]=2; cooColIndexHostPtr[4]=0; cooValHostPtr[4]=5.0; 
cooRowIndexHostPtr[5]=2; cooColIndexHostPtr[5]=2; cooValHostPtr[5]=6.0; 
cooRowIndexHostPtr[6]=2; cooColIndexHostPtr[6]=3; cooValHostPtr[6]=7.0; 
cooRowIndexHostPtr[7]=3; cooColIndexHostPtr[7]=1; cooValHostPtr[7]=8.0; 
cooRowIndexHostPtr[8]=3; cooColIndexHostPtr[8]=3; cooValHostPtr[8]=9.0;

Выведем матрицу на экран:

printf("Input data:\n"); 
for (i=0; i<nnz; i++){ 
printf("cooRowIndexHostPtr[%d]=%d ",i,cooRowIndexHostPtr[i]); 
printf("cooColIndexHostPtr[%d]=%d ",i,cooColIndexHostPtr[i]); 
printf("cooValHostPtr[%d]=%f \n",i,cooValHostPtr[i]); 
}

Создадим разреженный и плотный векторы:

// xVal= [100.0 200.0 400.0] (разреженный) 
// xInd= [0 1 3] 
// y = [10.0 20.0 30.0 40.0 | 50.0 60.0 70.0 80.0] (плотный) 
nnz_vector = 3; 
xIndHostPtr = (int *) malloc(nnz_vector*sizeof(xIndHostPtr[0])); 
xValHostPtr = (double *)malloc(nnz_vector*sizeof(xValHostPtr[0])); 
yHostPtr = (double *)malloc(2*n *sizeof(yHostPtr[0])); 
zHostPtr = (double *)malloc(2*(n+1) *sizeof(zHostPtr[0])); 
if((!xIndHostPtr) || (!xValHostPtr) || (!yHostPtr) || (!zHostPtr)){ 
CLEANUP("Host malloc failed (vectors)"); 
return EXIT_FAILURE; 
} 
yHostPtr[0] = 10.0; xIndHostPtr[0]=0; xValHostPtr[0]=100.0; 
yHostPtr[1] = 20.0; xIndHostPtr[1]=1; xValHostPtr[1]=200.0; 
yHostPtr[2] = 30.0; 
yHostPtr[3] = 40.0; xIndHostPtr[2]=3; xValHostPtr[2]=400.0; 
yHostPtr[4] = 50.0; 
yHostPtr[5] = 60.0; 
yHostPtr[6] = 70.0; 
yHostPtr[7] = 80.0;

Выведем векторы на экран:

for (j=0; j<2; j++){ 
    for (i=0; i<n; i++){ 
        printf("yHostPtr[%d,%d]=%f\n",i,j,yHostPtr[i+n*j]); 
    } 
} 
for (i=0; i<nnz_vector; i++){ 
printf("xIndHostPtr[%d]=%d ",i,xIndHostPtr[i]); 
printf("xValHostPtr[%d]=%f\n",i,xValHostPtr[i]); 
}

Выделим память на GPU и скопируем матрицы и векторы в нее:

cudaStat1 = cudaMalloc((void**)&cooRowIndex,nnz*sizeof(cooRowIndex[0])); 
cudaStat2 = cudaMalloc((void**)&cooColIndex,nnz*sizeof(cooColIndex[0])); 
cudaStat3 = cudaMalloc((void**)&cooVal, nnz*sizeof(cooVal[0])); 
cudaStat4 = cudaMalloc((void**)&y, 2*n*sizeof(y[0])); 
cudaStat5 = cudaMalloc((void**)&xInd,nnz_vector*sizeof(xInd[0])); 
cudaStat6 = cudaMalloc((void**)&xVal,nnz_vector*sizeof(xVal[0])); 
if ((cudaStat1 != cudaSuccess) || 
(cudaStat2 != cudaSuccess) || 
(cudaStat3 != cudaSuccess) || 
(cudaStat4 != cudaSuccess) || 
(cudaStat5 != cudaSuccess) || 
(cudaStat6 != cudaSuccess)) { 
CLEANUP("Device malloc failed"); 
return EXIT_FAILURE; 
} 
cudaStat1 = cudaMemcpy(cooRowIndex, cooRowIndexHostPtr, (size_t)(nnz*sizeof(cooRowIndex[0])), cudaMemcpyHostToDevice); 
cudaStat2 = cudaMemcpy(cooColIndex, cooColIndexHostPtr, (size_t)(nnz*sizeof(cooColIndex[0])), cudaMemcpyHostToDevice); 
cudaStat3 = cudaMemcpy(cooVal, cooValHostPtr, (size_t)(nnz*sizeof(cooVal[0])), cudaMemcpyHostToDevice); 
cudaStat4 = cudaMemcpy(y, yHostPtr, (size_t)(2*n*sizeof(y[0])), cudaMemcpyHostToDevice); 
cudaStat5 = cudaMemcpy(xInd, xIndHostPtr, (size_t)(nnz_vector*sizeof(xInd[0])), cudaMemcpyHostToDevice); 
cudaStat6 = cudaMemcpy(xVal, xValHostPtr, (size_t)(nnz_vector*sizeof(xVal[0])), cudaMemcpyHostToDevice); 
if ((cudaStat1 != cudaSuccess) || 
(cudaStat2 != cudaSuccess) || 
(cudaStat3 != cudaSuccess) || 
(cudaStat4 != cudaSuccess) || 
(cudaStat5 != cudaSuccess) || 
(cudaStat6 != cudaSuccess)) { 
CLEANUP("Memcpy from Host to Device failed"); 
return EXIT_FAILURE; 
}

Инициализируем библиотеку CUSPARSE:

status= CUSPARSECreate(&handle); 
if (status != CUSPARSE_STATUS_SUCCESS) { 
CLEANUP("CUSPARSE Library initialization failed"); 
return EXIT_FAILURE; 
}

Создадим и настроим дескриптор матрицы:

status= CUSPARSECreateMatDescr(&descra); 
if (status != CUSPARSE_STATUS_SUCCESS) { 
CLEANUP("Matrix descriptor initialization failed"); 
return EXIT_FAILURE; 
} 
CUSPARSESetMatType(descra,CUSPARSE_MATRIX_TYPE_GENERAL); 
CUSPARSESetMatIndexBase(descra,CUSPARSE_INDEX_BASE_ZERO);

Процедура преобразования (преобразование матрицы из формата COO в формат CSR)?

cudaStat1 = cudaMalloc((void**)&csrRowPtr,(n+1)*sizeof(csrRowPtr[0])); 
if (cudaStat1 != cudaSuccess) { 
CLEANUP("Device malloc failed (csrRowPtr)"); 
return EXIT_FAILURE; 
} 
status= CUSPARSEXcoo2csr(handle,cooRowIndex,nnz,n, 
csrRowPtr,CUSPARSE_INDEX_BASE_ZERO); 
if (status != CUSPARSE_STATUS_SUCCESS) { 
CLEANUP("Conversion from COO to CSR format failed"); 
return EXIT_FAILURE; 
} 
//csrRowPtr = [0 3 4 7 9]

Умножим разреженный вектор на плотный:

status= CUSPARSEDsctr(handle, nnz_vector, xVal, xInd, &y[n], CUSPARSE_INDEX_BASE_ZERO); 
if (status != CUSPARSE_STATUS_SUCCESS) { 
CLEANUP("Scatter from sparse to dense vector failed"); 
return EXIT_FAILURE; 
} 
//y = [10 20 30 40 | 100 200 70 400]

Умножим разреженную матрицу на плотный вектор:

status= CUSPARSEDcsrmv(handle,CUSPARSE_OPERATION_NON_TRANSPOSE, n, n, 2.0, descra, cooVal, csrRowPtr, cooColIndex, &y[0], 3.0, &y[n]); 
if (status != CUSPARSE_STATUS_SUCCESS) { 
CLEANUP("Matrix-vector multiplication failed"); 
return EXIT_FAILURE; 
}

Выведем на экран промежуточные результаты (у):

//y = [10 20 30 40 | 680 760 1230 2240] 
cudaMemcpy(yHostPtr, y, (size_t)(2*n*sizeof(y[0])), cudaMemcpyDeviceToHost); 
printf("Intermediate results:\n"); 
for (j=0; j<2; j++){ 
    for (i=0; i<n; i++){ 
    printf("yHostPtr[%d,%d]=%f\n",i,j,yHostPtr[i+n*j]); 
    } 
}

Умножим разреженную матрицу на плотный вектор-столбец):

cudaStat1 = cudaMalloc((void**)&z, 2*(n+1)*sizeof(z[0])); 
if (cudaStat1 != cudaSuccess) { 
CLEANUP("Device malloc failed (z)"); 
return EXIT_FAILURE; 
} 
cudaStat1 = cudaMemset((void *)z,0, 2*(n+1)*sizeof(z[0])); 
if (cudaStat1 != cudaSuccess) { 
CLEANUP("Memset on Device failed"); 
return EXIT_FAILURE; 
} 
status= CUSPARSEDcsrmm(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, n, 2, n, 5.0, descra, cooVal, csrRowPtr, cooColIndex, y, n, 0.0, z, n+1); 
if (status != CUSPARSE_STATUS_SUCCESS) { 
CLEANUP("Matrix-matrix multiplication failed"); 
return EXIT_FAILURE; 
}

Выведем конечные результаты на экран (z):

cudaStat1 = cudaMemcpy(zHostPtr, z, (size_t)(2*(n+1)*sizeof(z[0])), cudaMemcpyDeviceToHost); 
if (cudaStat1 != cudaSuccess) { 
CLEANUP("Memcpy from Device to Host failed"); 
return EXIT_FAILURE; 
} 
//z = [950 400 2550 2600 0 | 49300 15200 132300 131200 0] 
printf("Final results:\n"); 
for (j=0; j<2; j++){ 
for (i=0; i<n+1; i++){ 
printf("z[%d,%d]=%f\n",i,j,zHostPtr[i+(n+1)*j]); 
} 
}

Проверка результатов. Помните что CLEANUP() содержит вызов CUSPARSEDestroy(handle):

if ((zHostPtr[0] != 950.0) || 
(zHostPtr[1] != 400.0) || 
(zHostPtr[2] != 2550.0) || 
(zHostPtr[3] != 2600.0) || 
(zHostPtr[4] != 0.0) || 
(zHostPtr[5] != 49300.0) || 
(zHostPtr[6] != 15200.0) || 
(zHostPtr[7] != 132300.0) || 
(zHostPtr[8] != 131200.0) || 
(zHostPtr[9] != 0.0) || 
(yHostPtr[0] != 10.0) || 
(yHostPtr[1] != 20.0) || 
(yHostPtr[2] != 30.0) || 
(yHostPtr[3] != 40.0) || 
(yHostPtr[4] != 680.0) || 
(yHostPtr[5] != 760.0) || 
(yHostPtr[6] != 1230.0) || 
(yHostPtr[7] != 2240.0)){ 
CLEANUP("example test FAILED"); 
return EXIT_FAILURE; 
} 
else{ 
CLEANUP("example test PASSED"); 
return EXIT_SUCCESS; 
} 
}