我想要评估一个标量函数f(x),其中x是一个k维矢量(即f:R^k-> R)。在评估过程中,我必须执行许多矩阵运算:反演,乘法和寻找矩阵行列式和中等大小矩阵(大部分小于30x30)的轨迹。现在我想通过在GPU上使用不同的线程来同时在许多不同的xs上评估函数。这就是为什么我需要设备API。用cublas设备API计算矩阵行列式
我已经写了下面的代码来测试由cublas设备API cublasSgetrfBatched计算矩阵行列式,其中我首先找到矩阵的LU分解并计算U矩阵中所有对角线元素的乘积。我已经使用cublas返回的结果在GPU线程和CPU上完成了此操作。但是,GPU的结果没有任何意义,而CPU的结果是正确的。我用cuda-memcheck,但没有发现错误。有人可以帮助解决这个问题吗?非常感谢。
cat test2.cu
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
__host__ __device__ unsigned int IDX(unsigned int i,unsigned int j,unsigned int ld){return j*ld+i;}
#define PERR(call) \
if (call) {\
fprintf(stderr, "%s:%d Error [%s] on "#call"\n", __FILE__, __LINE__,\
cudaGetErrorString(cudaGetLastError()));\
exit(1);\
}
#define ERRCHECK \
if (cudaPeekAtLastError()) { \
fprintf(stderr, "%s:%d Error [%s]\n", __FILE__, __LINE__,\
cudaGetErrorString(cudaGetLastError()));\
exit(1);\
}
__device__ float
det_kernel(float *a_copy,unsigned int *n,cublasHandle_t *hdl){
int *info = (int *)malloc(sizeof(int));info[0]=0;
int batch=1;int *p = (int *)malloc(*n*sizeof(int));
float **a = (float **)malloc(sizeof(float *));
*a = a_copy;
cublasStatus_t status=cublasSgetrfBatched(*hdl, *n, a, *n, p, info, batch);
unsigned int i1;
float res=1;
for(i1=0;i1<(*n);++i1)res*=a_copy[IDX(i1,i1,*n)];
return res;
}
__global__ void runtest(float *a_i,unsigned int n){
cublasHandle_t hdl;cublasCreate_v2(&hdl);
printf("det on GPU:%f\n",det_kernel(a_i,&n,&hdl));
cublasDestroy_v2(hdl);
}
int
main(int argc, char **argv)
{
float a[] = {
1, 2, 3,
0, 4, 5,
1, 0, 0};
cudaSetDevice(1);//GTX780Ti on my machine,0 for GTX1080
unsigned int n=3,nn=n*n;
printf("a is \n");
for (int i = 0; i < n; ++i){
for (int j = 0; j < n; j++) printf("%f, ",a[IDX(i,j,n)]);
printf("\n");}
float *a_d;
PERR(cudaMalloc((void **)&a_d, nn*sizeof(float)));
PERR(cudaMemcpy(a_d, a, nn*sizeof(float), cudaMemcpyHostToDevice));
runtest<<<1, 1>>>(a_d,n);
cudaDeviceSynchronize();
ERRCHECK;
PERR(cudaMemcpy(a, a_d, nn*sizeof(float), cudaMemcpyDeviceToHost));
float res=1;
for (int i = 0; i < n; ++i)res*=a[IDX(i,i,n)];
printf("det on CPU:%f\n",res);
}
nvcc -arch=sm_35 -rdc=true -o test test2.cu -lcublas_device -lcudadevrt
./test
a is
1.000000, 0.000000, 1.000000,
2.000000, 4.000000, 0.000000,
3.000000, 5.000000, 0.000000,
det on GPU:0.000000
det on CPU:-2.000000
非常感谢,罗伯特。您的建议有效并可以产生预期结果。我有另一个设备例程,首先使用cublasSgetrfBatched获取LU分解,然后cublasSgetriBatched从LU输出获取逆。在他们之间,我需要使用cudaDeviceSynchronize()吗?对于小矩阵(3乘3),结果似乎相同。 – Jack
另一个问题是关于释放malloc在设备例程中创建的内存。如果我使用free()来释放内存,那么当代码使用-lcublas_device编译时(如果没有它),cuda-memcheck会产生错误。你碰巧有什么想法来解决这个问题吗?如果我没有释放记忆,会有什么后果?多谢。 – Jack
关于第一个问题,cuda [stream semantics](http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#streams)也适用于设备端操作。没有指定流的设备端启动应该启动到默认流(每个线程!),这意味着一个设备cublas操作A跟一个由同一个线程发出的设备cublas操作B应该被串行化。 B应该在A完成之后才开始。关于第二个问题,我需要确切知道你在做什么'free()'操作。 –