一、稀疏矩阵稀疏矩阵如果对线代不太了解的可能得回去看看相关的书籍。什么样的矩阵是稀疏矩阵呢其实很好理解就是字面意思。矩阵中稀稀拉拉的有几个非0的元素看上去就像地里没有长几棵庄稼。不过与之相对应的还有一个稠密矩阵就是矩阵中非0的数据太多了。很多人可能接触过数据压缩如音视频等其实就是把相同的元素用某种方式表示出来而不是简单的一一的用数据进行映射。稀疏矩阵也是如此一个特别大的稀疏矩阵可能需要一个巨大的数组来存储但其中真正有用的数据没有多少大量的空间其实被浪费。所以其除了可以进行原始方式存储外也可以使用各种压缩方式进行存储。包括常见的压缩稀疏行CSR格式和压缩稀疏列CSC格式。二、扫描的应用稀疏矩阵在实际的应用中非常常见如社交类软件的互相联系人AI中的机器学习和自然语言处理等中都存在这种应用。在前面学习过CUDA中的扫描和归约。其中就提到了对矩阵的乘法和加法的处理。所以说稀疏矩阵是对扫描应用的一个非常好的例子。特别是在迭代处理相关操作时更能体现出CUDA中的扫描的应用特点。CUDA中扫描有两种主要方式即Belloch算法和Hillis-steele算法不清楚的可以翻看前面的文章。其中Hillis-steele算法在step complexity方面占有优势而Belloch算法在work complexity有较大的优势。但前者更适合于中小规模的数据处理。所以对于类似稀疏矩阵这种处理一般使用Belloch算法。SpMVSparse Matrix-Vector Multiplication即稀疏矩阵向量乘法。它是科学计算、工程模拟以及大规模数据处理中应用最广泛、最基础的运算之一。三、例程下面看一个SpMV的例程#includecuda_runtime.h#includedevice_launch_parameters.h#includestdio.h#includeiostream#includestdlib.h#includemath.htypedeffloatreal;__device__ realreduceSum(real val,intsize){#pragmaunrollfor(intoffset16;offset0;offset1){val__shfl_down_sync(0xffffffff,val,offset,size);}returnval;}__global__voidspmvPerRow(intm,constint*rowPtr,constint*colIdx,constreal*values,constreal*x,real*y){intwarp_id(blockIdx.x*blockDim.xthreadIdx.x)/32;intlanethreadIdx.x31;introwwarp_id;if(rowm){introwStartrowPtr[row];introwEndrowPtr[row1];intnnzRowrowEnd-rowStart;real prod0.0f;if(lanennzRow){intidxrowStartlane;prodvalues[idx]*x[colIdx[idx]];}real sumreduceSum(prod,nnzRow);if(lane0){y[row]sum;}}}intmain(){constintm10000;constintn10000;constintnnzPerRow10;constintnnzm*nnzPerRow;// 100000int*hRowPtr(int*)malloc((m1)*sizeof(int));int*hColIdx(int*)malloc(nnz*sizeof(int));real*hValues(real*)malloc(nnz*sizeof(real));real*hX(real*)malloc(n*sizeof(real));real*hYCpu(real*)malloc(m*sizeof(real));real*hYGpu(real*)malloc(m*sizeof(real));if(!hRowPtr||!hColIdx||!hValues||!hX||!hYCpu||!hYGpu){std::cerrhost malloc err!std::endl;return-1;}hRowPtr[0]0;for(inti0;im;i){hRowPtr[i1]hRowPtr[i]nnzPerRow;}intused[10000]{0};for(inti0;im;i){for(intk0;kn;k){used[k]0;}introwStarthRowPtr[i];for(intk0;knnzPerRow;k){intcol;do{colrand()%n;}while(used[col]);used[col]1;hColIdx[rowStartk]col;hValues[rowStartk](real)rand()/RAND_MAX;}for(inti0;in;i){hX[i](real)rand()/RAND_MAX;}}int*dRowPtr,*dColIdx;real*dValues,*dX,*dY;cudaMalloc(dRowPtr,(m1)*sizeof(int));cudaMalloc(dColIdx,nnz*sizeof(int));cudaMalloc(dValues,nnz*sizeof(real));cudaMalloc(dX,n*sizeof(real));cudaMalloc(dY,m*sizeof(real));cudaMemcpy(dRowPtr,hRowPtr,(m1)*sizeof(int),cudaMemcpyHostToDevice);cudaMemcpy(dColIdx,hColIdx,nnz*sizeof(int),cudaMemcpyHostToDevice);cudaMemcpy(dValues,hValues,nnz*sizeof(real),cudaMemcpyHostToDevice);cudaMemcpy(dX,hX,n*sizeof(real),cudaMemcpyHostToDevice);constintthreadsPerBlock128;constintwarpsPerBlockthreadsPerBlock/32;constintblocksPerGrid(mwarpsPerBlock-1)/warpsPerBlock;spmvPerRowblocksPerGrid,threadsPerBlock(m,dRowPtr,dColIdx,dValues,dX,dY);cudaError_terrcudaGetLastError();if(err!cudaSuccess){std::cerrKernel launch failed: %scudaGetErrorString(err)std::endl;return-1;}free(hRowPtr);free(hColIdx);free(hValues);free(hX);free(hYCpu);free(hYGpu);cudaFree(dRowPtr);cudaFree(dColIdx);cudaFree(dValues);cudaFree(dX);cudaFree(dY);return0;}在上面的代码中设置了一个10000*10000的稀疏矩阵用Blloch算法进行处理。在前面的分析(并行任务和Warp分析)中提到过在共享内存的处理中一个重要的提醒就是要防范出现Bank conflict。这就需要开发者使用一些数据处理的技巧来进行控制对内存的访问了。四、总结矩阵处理是大数据、科学计算以及现在大火的AI的基础数据运算的重要基础之一。想搞AI离不开对矩阵的相关处理的学习和分析能力。虽然现在大多数的语言和框架中都提供了对其的封装但如果没有掌握其原理和内容的话应用起来还是存在着卡点的。这也是前面为什么在路线里提出要掌握数学知识的原因。