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

Question

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

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

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