c++ - 使用 cublas 设备 API 计算矩阵行列式

Question

我正在尝试计算标量函数 f(x)，其中 x 是 k 维 vector (即 f:R^k->R)。在评估期间，我必须执行许多矩阵运算:求逆、乘法和查找中等大小矩阵(大多数小于 30x30)的矩阵行列式和迹。现在我想通过在 GPU 上使用不同的线程同时在许多不同的 xs 上评估函数。这就是我需要设备 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

Answer 1

cublas 设备调用是异步的。

这意味着它们在 cublas 调用完成之前将控制权返回给调用线程。

如果您希望调用线程能够直接处理结果(就像您在此处计算 res 时所做的那样)，则必须在开始计算之前强制同步以等待结果。

在主机端计算中看不到这一点，因为在父内核终止之前，任何设备事件(包括 cublas 设备动态并行性)都存在隐式同步。

因此，如果您在设备 cublas 调用之后添加同步，如下所示:

cublasStatus_t status=cublasSgetrfBatched(*hdl, *n, a, *n, p, info, batch); 
cudaDeviceSynchronize(); // add this line

我想您会看到设备计算和主机计算之间的匹配，正如您所期望的那样。