与 cusparse 相比，cublas 异常缓慢
cublas is unusually slow compare to cusparse
我正在尝试 运行 一些测试来比较不同稀疏度（使用 Titan X）下的 cusparse 和 cublas 性能，这里是名为 "testcusparsevector.cpp" 的主要代码：
#include <stdio.h>
#include <iostream>
#include <vector>
#include <cstdlib>
#include <fstream>
#include <time.h>
#include <cuda_runtime.h>
#include <cublas.h>
#include <cusparse_v2.h>
#include <cublas_v2.h>
#include <assert.h>
#define M 6
#define N 5
#define IDX2C(i,j,ld) (((j)*(ld))+(i))
// /home/gpu1/Install/OpenBLAS-0.2.14
#define CHECK_EQ(a,b) do { \
if ((a) != (b)) { \
cout <<__FILE__<<" : "<< __LINE__<<" : check failed because "<<a<<"!="<<b<<endl;\
exit(1);\
}\
} while(0)
#define CUBLAS_CHECK(condition) \
do {\
cublasStatus_t status = condition; \
CHECK_EQ(status, CUBLAS_STATUS_SUCCESS); \
} while(0)
#define CUSPARSE_CHECK(condition)\
do {\
cusparseStatus_t status = condition; \
switch(status)\
{\
case CUSPARSE_STATUS_NOT_INITIALIZED:\
cout << "CUSPARSE_STATUS_NOT_INITIALIZED" << endl;\
break;\
case CUSPARSE_STATUS_ALLOC_FAILED:\
cout << "CUSPARSE_STATUS_ALLOC_FAILED" << endl;\
break;\
case CUSPARSE_STATUS_INVALID_VALUE:\
cout << "CUSPARSE_STATUS_INVALID_VALUE" << endl;\
break;\
case CUSPARSE_STATUS_ARCH_MISMATCH:\
cout << "CUSPARSE_STATUS_ARCH_MISMATCH" << endl;\
break;\
case CUSPARSE_STATUS_MAPPING_ERROR:\
cout << "CUSPARSE_STATUS_MAPPING_ERROR" << endl;\
break;\
case CUSPARSE_STATUS_EXECUTION_FAILED:\
cout << "CUSPARSE_STATUS_EXECUTION_FAILED" << endl;\
break;\
case CUSPARSE_STATUS_INTERNAL_ERROR:\
cout << "CUSPARSE_STATUS_INTERNAL_ERROR" << endl;\
break;\
case CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED:\
cout << "CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED" << endl;\
break;\
case CUSPARSE_STATUS_ZERO_PIVOT:\
cout << "CUSPARSE_STATUS_ZERO_PIVOT" << endl;\
}\
CHECK_EQ(status, CUSPARSE_STATUS_SUCCESS); \
} while(0)
#define CUDA_CHECK(condition)\
do {\
cudaError_t error = condition;\
CHECK_EQ(error, cudaSuccess);\
} while(0)
//check after kernel function
#define CUDA_POST_KERNEL_CHECK CUDA_CHECK(cudaPeekAtLastError())
#define __TIMING__ 1
#if __TIMING__
#define INIT_TIMER cudaEvent_t start, stop; \
float milliseconds = 0; \
float sum = 0;\
cudaEventCreate( &start );\
cudaEventCreate( &stop );
#define TIC { cudaEventRecord( start ); }
#if __CUDNN__
#define PREDEFNAME "CUDNN"
#else
#define PREDEFNAME "CUDA"
#endif
#define TOC(a) { cudaEventRecord( stop ); \
cudaEventSynchronize( stop ); \
cudaEventElapsedTime( &milliseconds, start, stop ); \
printf( "GPU Execution time of %s_%s: %f ms\n",PREDEFNAME, a, milliseconds ); \
sum += milliseconds;\
fflush(stdout); }
#define CLOSE_TIMER {cudaEventDestroy(start); cudaEventDestroy(stop); }
#endif
using namespace std;
void dispArray(double* array, size_t width, size_t height) {
for (int i=0; i < height;i++ ) {
for (int j=0;j < width;j++) {
cout << array[j*height+i] << ' ';
}
cout << endl;
}
cout << endl;
}
int main()
{
srand(time(NULL));
const int num_loop = 1;
const int inside_loop = 1000;
// const int WIDTH = 512*3*3;
// const int HEIGHT = 512;
// const int WIDTHOUT = 36;
const int WIDTH = 4608;
const int HEIGHT = 512;
const int WIDTHOUT = 144;
// const int WIDTH = 18500;
// const int HEIGHT = 512;
// const int WIDTHOUT = 1;
// const int WIDTH = 3;
// const int HEIGHT = 5;
// const int WIDTHOUT = 2;
INIT_TIMER
ofstream myfile;
myfile.open("test_sparsity.log");
cudaError_t cudaStat;
cusparseStatus_t stat;
cusparseHandle_t handle;
cublasHandle_t handleblas;
double *devPtrOutput;
double *devPtrOutput2;
double *devPtrRand;
double *devPtrSec;
CUDA_CHECK(cudaMalloc((void **)&(devPtrOutput), sizeof(double)*HEIGHT*WIDTHOUT));
CUDA_CHECK(cudaMalloc((void **)&(devPtrOutput2), sizeof(double)*HEIGHT*WIDTHOUT));
CUDA_CHECK(cudaMalloc((void **)&(devPtrRand), sizeof(double)*WIDTH*WIDTHOUT));
CUDA_CHECK(cudaMalloc((void **)&(devPtrSec), sizeof(double)*WIDTH*HEIGHT));
const double alpha=1.0;
const double beta=0.0;
double *csrVal;
int *csrRowPtr;
int *csrColInd;
const bool SPARSE = true;
long a = clock();
long temp = clock();
cusparseMatDescr_t descr;
CUSPARSE_CHECK(cusparseCreateMatDescr(&descr));
cusparseSetMatType(descr, CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(descr, CUSPARSE_INDEX_BASE_ZERO);
int nnz;
CUSPARSE_CHECK(cusparseCreate(&handle));
CUBLAS_CHECK(cublasCreate(&handleblas));
int *nnzPerRow_gpu;
CUDA_CHECK(cudaMalloc((void **)&(nnzPerRow_gpu), sizeof(int)*HEIGHT));
CUDA_CHECK(cudaMalloc((void **)&(csrRowPtr), sizeof(int)*(HEIGHT+1)));
double density_array[1] = {0.9999};//, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1 ,0.09, 0.08, 0.07, 0.06, 0.05 ,0.04, 0.03, 0.02, 0.01};
for (int inddense=0;inddense < 1;inddense++) {
double DENSITY = density_array[inddense];
int num_non_zeros = DENSITY * (WIDTH * HEIGHT);
CUDA_CHECK(cudaMalloc((void **)&(csrColInd), sizeof(int)*num_non_zeros));
CUDA_CHECK(cudaMalloc((void **)&(csrVal), sizeof(double)*num_non_zeros));
INIT_TIMER
for (int iter=0; iter < num_loop;iter++) {
vector<double> randVec(WIDTH*WIDTHOUT, 0);
vector<double> secArray(WIDTH*HEIGHT, 0);
vector<int> temp(WIDTH*HEIGHT, 1);
for (int j = 0; j < WIDTH*WIDTHOUT; j++) {
randVec[j]=(double)(rand()%100000)/100;
}
for (int x, i = 0; i < num_non_zeros;i++) {
do
{
x = rand() % (WIDTH*HEIGHT);
} while(temp[x] == 0);
temp[x]=0;
secArray[x]=(double)(rand()%100000)/100;
}
int count = 0;
for(int i=0;i < WIDTH*HEIGHT;i++) {
if (secArray[i] != 0) {
count++;
}
}
// randVec = {2,2,2,3,3,3};
// secArray = {0,5,0,2,5,8,7,0,0,0,0,2,0,4,4};
CUDA_CHECK(cudaMemcpy(devPtrRand, &randVec[0], sizeof(double)*WIDTH*WIDTHOUT, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(devPtrSec, &secArray[0], sizeof(double)*WIDTH*HEIGHT, cudaMemcpyHostToDevice));
if (SPARSE) {
CUSPARSE_CHECK(cusparseDnnz(handle, CUSPARSE_DIRECTION_ROW, HEIGHT, WIDTH, descr, devPtrSec, HEIGHT, nnzPerRow_gpu, &nnz));
CUSPARSE_CHECK(cusparseDdense2csr(handle, HEIGHT, WIDTH, descr,devPtrSec,HEIGHT,nnzPerRow_gpu,csrVal,csrRowPtr,csrColInd));
}
// vector<double> tempcsrVal(nnz,0);
// vector<int> tempcsrRowPtr(HEIGHT+1);
// vector<int> tempcsrColInd(nnz,0);
// CUDA_CHECK(cudaMemcpy(&tempcsrVal[0], csrVal, sizeof(double)*nnz, cudaMemcpyDeviceToHost));
// CUDA_CHECK(cudaMemcpy(&tempcsrRowPtr[0], csrRowPtr, sizeof(int)*(HEIGHT+1), cudaMemcpyDeviceToHost));
// CUDA_CHECK(cudaMemcpy(&tempcsrColInd[0], csrColInd, sizeof(int)*nnz, cudaMemcpyDeviceToHost));
// for (int i =0; i < nnz;i++) {
// cout << tempcsrVal[i] << " ";
// }
// cout << endl;
// for (int i =0; i < HEIGHT+1;i++) {
// cout << tempcsrRowPtr[i] << " ";
// }
// cout << endl;
// for (int i =0; i < nnz;i++) {
// cout << tempcsrColInd[i] << " ";
// }
// cout << endl;
cudaDeviceSynchronize();
TIC
for (int i=0 ; i < inside_loop;i++) {
if (WIDTHOUT == 1) {
// TIC
CUSPARSE_CHECK(cusparseDcsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
HEIGHT, WIDTH, nnz, &alpha, descr, csrVal, csrRowPtr, csrColInd,
devPtrRand, &beta, devPtrOutput));
// TOC("csrmv")
} else {
// TIC
CUSPARSE_CHECK(cusparseDcsrmm(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
HEIGHT, WIDTHOUT, WIDTH, nnz, &alpha, descr, csrVal, csrRowPtr,
csrColInd, devPtrRand, WIDTH, &beta, devPtrOutput, HEIGHT));
// TOC("csrmm")
}
}
TOC("csr")
TIC
for (int i=0 ; i < inside_loop;i++) {
if (WIDTHOUT == 1) {
// TIC
CUBLAS_CHECK(cublasDgemv(handleblas, CUBLAS_OP_N, HEIGHT, WIDTH, &alpha, devPtrSec, HEIGHT , devPtrRand, 1, &beta, devPtrOutput2, 1));
// TOC("dgemv")
} else {
// TIC
CUBLAS_CHECK(cublasDgemm(handleblas, CUBLAS_OP_N, CUBLAS_OP_N, HEIGHT, WIDTHOUT, WIDTH, &alpha, devPtrSec, HEIGHT, devPtrRand, WIDTH, &beta, devPtrOutput2, HEIGHT));
// TOC("dgemm")
}
}
TOC("blas")
#if 0
vector<double> output(HEIGHT*WIDTHOUT, 0);
vector<double> output2(HEIGHT*WIDTHOUT, 0);
CUDA_CHECK(cudaMemcpy(&output[0], devPtrOutput, sizeof(double)*HEIGHT*WIDTHOUT, cudaMemcpyDeviceToHost));
CUDA_CHECK(cudaMemcpy(&output2[0], devPtrOutput2, sizeof(double)*HEIGHT*WIDTHOUT, cudaMemcpyDeviceToHost));
dispArray(&output[0], WIDTHOUT, HEIGHT);
cout << endl;
for (int i=0;i < WIDTHOUT * HEIGHT;i++) {
if (output[i] != output2[i]) {
cout << "error: " << i << " " << (output[i] - output2[i]) << " " << output[i] << endl;
}
}
#endif
}
cout << DENSITY << " " << sum/num_loop << endl;
myfile << DENSITY << " " << sum/num_loop << endl;
cudaFree(csrColInd);
cudaFree(csrVal);
}
myfile.close();
cudaFree(csrRowPtr);
cudaFree(devPtrOutput);
cudaFree(devPtrRand);
cudaFree(devPtrSec);
}
但是在使用
编译代码之后
g++ -std=c++1y -O3 -I/usr/local/cuda/include -o testcusparsevector testcusparsevector.cpp -L/usr/local/cuda/lib64 -lcudart -lcublas -lcusparse
这是输出：
GPU Execution time of CUDA_csr: 4818.447266 ms
GPU Execution time of CUDA_blas: 5024.459961 ms
这应该意味着即使我的密度是 0.999，cusparseDcsrmm 仍然比 cublasDgemm 快，我已经检查了好的结果，并且与其他示例相比，似乎问题来自 cublas 这是太慢了。
你知道它是从哪里来的吗？
编辑：我试图将值更改为浮点数，结果更符合我的要求，显然，cublas 不是为双重计算而设计的...
提前致谢。
Titan X（以及 maxwell GPU 系列的所有当前成员）在双精度浮点运算和单精度浮点运算之间的吞吐量比率为 1:32。
通常，稀疏矩阵运算受内存带宽限制，而密集矩阵-矩阵乘法是计算受限问题的一个示例。
因此，在您的示例中，您正在处理一个通常受计算限制的问题，并且 运行 它作为一个稀疏矩阵乘法在具有相对大量内存带宽和相对少量双精度计算吞吐量。
这种情况会导致两个 API 之间的界限变得模糊，而 CUBLAS API 通常会更快地进行这种比较。
如果您将代码切换为使用 float 而不是我认为您已经尝试过的 double，您将看到 CUBLAS 再次获胜。同样，如果您 运行 在单精度和双精度吞吐量之间具有不同比率的 GPU 上按原样编写代码，您也会看到 CUBLAS 在那里再次获胜。
apparently, cublas is not made for double computation...
与其这样说，不如说 GTX Titan X 不是（主要）用于双重计算的。尝试使用 Tesla K80、K40 或其他具有更接近双吞吐量与单吞吐量比率的 GPU。
这是你的程序 运行 在 "unboosted" Tesla K40 上的输出：
$ ./testcusparsevector
GPU Execution time of CUDA_csr: 8870.386719 ms
GPU Execution time of CUDA_blas: 1045.211792 ms
免责声明：我没有尝试研究您的代码。我仔细看了看，没有发现任何明显的问题。但可能有些问题我没有发现。
我正在尝试 运行 一些测试来比较不同稀疏度（使用 Titan X）下的 cusparse 和 cublas 性能，这里是名为 "testcusparsevector.cpp" 的主要代码：
#include <stdio.h>
#include <iostream>
#include <vector>
#include <cstdlib>
#include <fstream>
#include <time.h>
#include <cuda_runtime.h>
#include <cublas.h>
#include <cusparse_v2.h>
#include <cublas_v2.h>
#include <assert.h>
#define M 6
#define N 5
#define IDX2C(i,j,ld) (((j)*(ld))+(i))
// /home/gpu1/Install/OpenBLAS-0.2.14
#define CHECK_EQ(a,b) do { \
if ((a) != (b)) { \
cout <<__FILE__<<" : "<< __LINE__<<" : check failed because "<<a<<"!="<<b<<endl;\
exit(1);\
}\
} while(0)
#define CUBLAS_CHECK(condition) \
do {\
cublasStatus_t status = condition; \
CHECK_EQ(status, CUBLAS_STATUS_SUCCESS); \
} while(0)
#define CUSPARSE_CHECK(condition)\
do {\
cusparseStatus_t status = condition; \
switch(status)\
{\
case CUSPARSE_STATUS_NOT_INITIALIZED:\
cout << "CUSPARSE_STATUS_NOT_INITIALIZED" << endl;\
break;\
case CUSPARSE_STATUS_ALLOC_FAILED:\
cout << "CUSPARSE_STATUS_ALLOC_FAILED" << endl;\
break;\
case CUSPARSE_STATUS_INVALID_VALUE:\
cout << "CUSPARSE_STATUS_INVALID_VALUE" << endl;\
break;\
case CUSPARSE_STATUS_ARCH_MISMATCH:\
cout << "CUSPARSE_STATUS_ARCH_MISMATCH" << endl;\
break;\
case CUSPARSE_STATUS_MAPPING_ERROR:\
cout << "CUSPARSE_STATUS_MAPPING_ERROR" << endl;\
break;\
case CUSPARSE_STATUS_EXECUTION_FAILED:\
cout << "CUSPARSE_STATUS_EXECUTION_FAILED" << endl;\
break;\
case CUSPARSE_STATUS_INTERNAL_ERROR:\
cout << "CUSPARSE_STATUS_INTERNAL_ERROR" << endl;\
break;\
case CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED:\
cout << "CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED" << endl;\
break;\
case CUSPARSE_STATUS_ZERO_PIVOT:\
cout << "CUSPARSE_STATUS_ZERO_PIVOT" << endl;\
}\
CHECK_EQ(status, CUSPARSE_STATUS_SUCCESS); \
} while(0)
#define CUDA_CHECK(condition)\
do {\
cudaError_t error = condition;\
CHECK_EQ(error, cudaSuccess);\
} while(0)
//check after kernel function
#define CUDA_POST_KERNEL_CHECK CUDA_CHECK(cudaPeekAtLastError())
#define __TIMING__ 1
#if __TIMING__
#define INIT_TIMER cudaEvent_t start, stop; \
float milliseconds = 0; \
float sum = 0;\
cudaEventCreate( &start );\
cudaEventCreate( &stop );
#define TIC { cudaEventRecord( start ); }
#if __CUDNN__
#define PREDEFNAME "CUDNN"
#else
#define PREDEFNAME "CUDA"
#endif
#define TOC(a) { cudaEventRecord( stop ); \
cudaEventSynchronize( stop ); \
cudaEventElapsedTime( &milliseconds, start, stop ); \
printf( "GPU Execution time of %s_%s: %f ms\n",PREDEFNAME, a, milliseconds ); \
sum += milliseconds;\
fflush(stdout); }
#define CLOSE_TIMER {cudaEventDestroy(start); cudaEventDestroy(stop); }
#endif
using namespace std;
void dispArray(double* array, size_t width, size_t height) {
for (int i=0; i < height;i++ ) {
for (int j=0;j < width;j++) {
cout << array[j*height+i] << ' ';
}
cout << endl;
}
cout << endl;
}
int main()
{
srand(time(NULL));
const int num_loop = 1;
const int inside_loop = 1000;
// const int WIDTH = 512*3*3;
// const int HEIGHT = 512;
// const int WIDTHOUT = 36;
const int WIDTH = 4608;
const int HEIGHT = 512;
const int WIDTHOUT = 144;
// const int WIDTH = 18500;
// const int HEIGHT = 512;
// const int WIDTHOUT = 1;
// const int WIDTH = 3;
// const int HEIGHT = 5;
// const int WIDTHOUT = 2;
INIT_TIMER
ofstream myfile;
myfile.open("test_sparsity.log");
cudaError_t cudaStat;
cusparseStatus_t stat;
cusparseHandle_t handle;
cublasHandle_t handleblas;
double *devPtrOutput;
double *devPtrOutput2;
double *devPtrRand;
double *devPtrSec;
CUDA_CHECK(cudaMalloc((void **)&(devPtrOutput), sizeof(double)*HEIGHT*WIDTHOUT));
CUDA_CHECK(cudaMalloc((void **)&(devPtrOutput2), sizeof(double)*HEIGHT*WIDTHOUT));
CUDA_CHECK(cudaMalloc((void **)&(devPtrRand), sizeof(double)*WIDTH*WIDTHOUT));
CUDA_CHECK(cudaMalloc((void **)&(devPtrSec), sizeof(double)*WIDTH*HEIGHT));
const double alpha=1.0;
const double beta=0.0;
double *csrVal;
int *csrRowPtr;
int *csrColInd;
const bool SPARSE = true;
long a = clock();
long temp = clock();
cusparseMatDescr_t descr;
CUSPARSE_CHECK(cusparseCreateMatDescr(&descr));
cusparseSetMatType(descr, CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(descr, CUSPARSE_INDEX_BASE_ZERO);
int nnz;
CUSPARSE_CHECK(cusparseCreate(&handle));
CUBLAS_CHECK(cublasCreate(&handleblas));
int *nnzPerRow_gpu;
CUDA_CHECK(cudaMalloc((void **)&(nnzPerRow_gpu), sizeof(int)*HEIGHT));
CUDA_CHECK(cudaMalloc((void **)&(csrRowPtr), sizeof(int)*(HEIGHT+1)));
double density_array[1] = {0.9999};//, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1 ,0.09, 0.08, 0.07, 0.06, 0.05 ,0.04, 0.03, 0.02, 0.01};
for (int inddense=0;inddense < 1;inddense++) {
double DENSITY = density_array[inddense];
int num_non_zeros = DENSITY * (WIDTH * HEIGHT);
CUDA_CHECK(cudaMalloc((void **)&(csrColInd), sizeof(int)*num_non_zeros));
CUDA_CHECK(cudaMalloc((void **)&(csrVal), sizeof(double)*num_non_zeros));
INIT_TIMER
for (int iter=0; iter < num_loop;iter++) {
vector<double> randVec(WIDTH*WIDTHOUT, 0);
vector<double> secArray(WIDTH*HEIGHT, 0);
vector<int> temp(WIDTH*HEIGHT, 1);
for (int j = 0; j < WIDTH*WIDTHOUT; j++) {
randVec[j]=(double)(rand()%100000)/100;
}
for (int x, i = 0; i < num_non_zeros;i++) {
do
{
x = rand() % (WIDTH*HEIGHT);
} while(temp[x] == 0);
temp[x]=0;
secArray[x]=(double)(rand()%100000)/100;
}
int count = 0;
for(int i=0;i < WIDTH*HEIGHT;i++) {
if (secArray[i] != 0) {
count++;
}
}
// randVec = {2,2,2,3,3,3};
// secArray = {0,5,0,2,5,8,7,0,0,0,0,2,0,4,4};
CUDA_CHECK(cudaMemcpy(devPtrRand, &randVec[0], sizeof(double)*WIDTH*WIDTHOUT, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(devPtrSec, &secArray[0], sizeof(double)*WIDTH*HEIGHT, cudaMemcpyHostToDevice));
if (SPARSE) {
CUSPARSE_CHECK(cusparseDnnz(handle, CUSPARSE_DIRECTION_ROW, HEIGHT, WIDTH, descr, devPtrSec, HEIGHT, nnzPerRow_gpu, &nnz));
CUSPARSE_CHECK(cusparseDdense2csr(handle, HEIGHT, WIDTH, descr,devPtrSec,HEIGHT,nnzPerRow_gpu,csrVal,csrRowPtr,csrColInd));
}
// vector<double> tempcsrVal(nnz,0);
// vector<int> tempcsrRowPtr(HEIGHT+1);
// vector<int> tempcsrColInd(nnz,0);
// CUDA_CHECK(cudaMemcpy(&tempcsrVal[0], csrVal, sizeof(double)*nnz, cudaMemcpyDeviceToHost));
// CUDA_CHECK(cudaMemcpy(&tempcsrRowPtr[0], csrRowPtr, sizeof(int)*(HEIGHT+1), cudaMemcpyDeviceToHost));
// CUDA_CHECK(cudaMemcpy(&tempcsrColInd[0], csrColInd, sizeof(int)*nnz, cudaMemcpyDeviceToHost));
// for (int i =0; i < nnz;i++) {
// cout << tempcsrVal[i] << " ";
// }
// cout << endl;
// for (int i =0; i < HEIGHT+1;i++) {
// cout << tempcsrRowPtr[i] << " ";
// }
// cout << endl;
// for (int i =0; i < nnz;i++) {
// cout << tempcsrColInd[i] << " ";
// }
// cout << endl;
cudaDeviceSynchronize();
TIC
for (int i=0 ; i < inside_loop;i++) {
if (WIDTHOUT == 1) {
// TIC
CUSPARSE_CHECK(cusparseDcsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
HEIGHT, WIDTH, nnz, &alpha, descr, csrVal, csrRowPtr, csrColInd,
devPtrRand, &beta, devPtrOutput));
// TOC("csrmv")
} else {
// TIC
CUSPARSE_CHECK(cusparseDcsrmm(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
HEIGHT, WIDTHOUT, WIDTH, nnz, &alpha, descr, csrVal, csrRowPtr,
csrColInd, devPtrRand, WIDTH, &beta, devPtrOutput, HEIGHT));
// TOC("csrmm")
}
}
TOC("csr")
TIC
for (int i=0 ; i < inside_loop;i++) {
if (WIDTHOUT == 1) {
// TIC
CUBLAS_CHECK(cublasDgemv(handleblas, CUBLAS_OP_N, HEIGHT, WIDTH, &alpha, devPtrSec, HEIGHT , devPtrRand, 1, &beta, devPtrOutput2, 1));
// TOC("dgemv")
} else {
// TIC
CUBLAS_CHECK(cublasDgemm(handleblas, CUBLAS_OP_N, CUBLAS_OP_N, HEIGHT, WIDTHOUT, WIDTH, &alpha, devPtrSec, HEIGHT, devPtrRand, WIDTH, &beta, devPtrOutput2, HEIGHT));
// TOC("dgemm")
}
}
TOC("blas")
#if 0
vector<double> output(HEIGHT*WIDTHOUT, 0);
vector<double> output2(HEIGHT*WIDTHOUT, 0);
CUDA_CHECK(cudaMemcpy(&output[0], devPtrOutput, sizeof(double)*HEIGHT*WIDTHOUT, cudaMemcpyDeviceToHost));
CUDA_CHECK(cudaMemcpy(&output2[0], devPtrOutput2, sizeof(double)*HEIGHT*WIDTHOUT, cudaMemcpyDeviceToHost));
dispArray(&output[0], WIDTHOUT, HEIGHT);
cout << endl;
for (int i=0;i < WIDTHOUT * HEIGHT;i++) {
if (output[i] != output2[i]) {
cout << "error: " << i << " " << (output[i] - output2[i]) << " " << output[i] << endl;
}
}
#endif
}
cout << DENSITY << " " << sum/num_loop << endl;
myfile << DENSITY << " " << sum/num_loop << endl;
cudaFree(csrColInd);
cudaFree(csrVal);
}
myfile.close();
cudaFree(csrRowPtr);
cudaFree(devPtrOutput);
cudaFree(devPtrRand);
cudaFree(devPtrSec);
}
但是在使用
编译代码之后g++ -std=c++1y -O3 -I/usr/local/cuda/include -o testcusparsevector testcusparsevector.cpp -L/usr/local/cuda/lib64 -lcudart -lcublas -lcusparse
这是输出：
GPU Execution time of CUDA_csr: 4818.447266 ms
GPU Execution time of CUDA_blas: 5024.459961 ms
这应该意味着即使我的密度是 0.999，cusparseDcsrmm 仍然比 cublasDgemm 快，我已经检查了好的结果，并且与其他示例相比，似乎问题来自 cublas 这是太慢了。
你知道它是从哪里来的吗？
编辑：我试图将值更改为浮点数，结果更符合我的要求，显然，cublas 不是为双重计算而设计的...
提前致谢。
Titan X（以及 maxwell GPU 系列的所有当前成员）在双精度浮点运算和单精度浮点运算之间的吞吐量比率为 1:32。
通常，稀疏矩阵运算受内存带宽限制，而密集矩阵-矩阵乘法是计算受限问题的一个示例。
因此，在您的示例中，您正在处理一个通常受计算限制的问题，并且 运行 它作为一个稀疏矩阵乘法在具有相对大量内存带宽和相对少量双精度计算吞吐量。
这种情况会导致两个 API 之间的界限变得模糊，而 CUBLAS API 通常会更快地进行这种比较。
如果您将代码切换为使用 float 而不是我认为您已经尝试过的 double，您将看到 CUBLAS 再次获胜。同样，如果您 运行 在单精度和双精度吞吐量之间具有不同比率的 GPU 上按原样编写代码，您也会看到 CUBLAS 在那里再次获胜。
apparently, cublas is not made for double computation...
与其这样说，不如说 GTX Titan X 不是（主要）用于双重计算的。尝试使用 Tesla K80、K40 或其他具有更接近双吞吐量与单吞吐量比率的 GPU。
这是你的程序 运行 在 "unboosted" Tesla K40 上的输出：
$ ./testcusparsevector
GPU Execution time of CUDA_csr: 8870.386719 ms
GPU Execution time of CUDA_blas: 1045.211792 ms
免责声明：我没有尝试研究您的代码。我仔细看了看，没有发现任何明显的问题。但可能有些问题我没有发现。