Рассчитать матричные детерминанты с помощью API устройства cublas

Question

Я пытаюсь оценить скалярную функцию f(x), где x является k-мерным вектором (т.е. f:R^k->R). Во время оценки мне нужно выполнить много матричных операций: инверсию, умножение и поиск матричных определителей и трасс для матриц среднего размера (большинство из них меньше 30x30). Теперь я хочу оценить функцию для множества разных х одновременно, используя разные потоки в графическом процессоре. Вот почему мне нужно устройство API.

Я написал следующие коды, чтобы проверить вычисление определителей матрицы с помощью API устройства cublas, cublasSgetrfBatched, где я сначала нахожу разложение матрицы LU и вычисляю произведение всех диагональных элементов в матрице U. Я сделал это как на GPU-потоке, так и на CPU, используя результат, возвращаемый cublas. Но результат от графического процессора не имеет никакого смысла, в то время как результат на процессоре правильный. Я использовал 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), прежде чем родительское ядро ​​завершит работу.

Так что, если вы добавите добавить синхронизацию после вызова устройства, как это:

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

Я думаю, вы увидите соответствие между вычислением устройства и вычислением хоста, как вы ожидаете.