用cublas設備API計算矩陣行列式

Question

我想要評估一個標量函數f（x），其中x是一個k維矢量（即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 

我想你會看到該設備的計算和主機計算之間的匹配，如您所願。