cublasデバイスAPIを使用して行列の行列式を計算する

Question

スカラー関数f（x）を評価しようとしています。ここで、xはk次元のベクトルです（つまり、f：R^k-> R）。評価の間、私は多くの行列演算を実行しなければなりません：逆行列、乗法、行列の行列式と中程度の大きさの行列のトレースを見つける必要があります（大部分は30x30未満です）。今私は、GPU上の異なるスレッドを使用して、同時に多くの異なるxの関数を評価したいと思います。だから私はデバイスのAPIが必要です。

私はcublasデバイスAPI、cublasSgetrfBatchedによって行列の行列式を計算するテストを行うために次のコードを書いています。まず、行列のLU分解を見つけ、U行列のすべての対角要素の積を計算します。私は、cublasによって返された結果を使って、GPUスレッドとCPUの両方でこれを行っています。しかし、CPU上の結果が正しい間は、GPUの結果は意味をなさない。私は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 

私はあなたが期待どおりには、デバイスの計算とホスト計算との間の一致を見ると思います。