matrice cuBLAS inverse beaucoup plus lent que MATLAB

Question

Dans mon projet actuel, j'essaie de calculer l'inverse d'une grande matrice (n> 2000) avec cuBLAS. Le calcul inverse est effectué, mais pour une raison quelconque, les temps de calcul sont significativement plus lents comparés à ceux effectués dans MATLAB.

J'ai joint un exemple de calcul effectué sur des matrices aléatoires en utilisant ma mise en œuvre dans les deux langues ainsi que les résultats de performance.

Toute aide ou suggestion sur ce qui pourrait causer ce ralentissement serait grandement appréciée.

Merci d'avance.

Comparaison

cuBLAS vs Matlab

N = 500: cuBLAS ~ 0,130 sec, ~ 0,066 sec Matlab -> ~ 1.97x plus lent

N = 1000: cuBLAS ~ 0.898 sec, MATLAB ~ 0.311 sec -> ~ 2.89x plus lent

N = 2000: cuBLAS ~ 6.667 sec, MATLAB ~ 0.659 sec -> ~ 10.12x plus lent

N = 4000: cuBLAS ~ 51,860 sec, MATLAB ~ 4,296 s -> ~ 12.07x lente

C++ code

#include <string> 
#include <cuda_runtime.h> 
#include <cublas_v2.h> 
#include <conio.h> 

#define CUDA_CALL(res, str) { if (res != cudaSuccess) { printf("CUDA Error : %s : %s %d : ERR %s\n", str, __FILE__, __LINE__, cudaGetErrorName(res)); } } 
#define CUBLAS_CALL(res, str) { if (res != CUBLAS_STATUS_SUCCESS) { printf("CUBLAS Error : %s : %s %d : ERR %d\n", str, __FILE__, __LINE__, int(res)); } } 

static cudaEvent_t cu_TimerStart; 
static cudaEvent_t cu_TimerStop; 

void d_CUDATimerStart(void) 
{ 
    CUDA_CALL(cudaEventCreate(&cu_TimerStart), "Failed to create start event!"); 
    CUDA_CALL(cudaEventCreate(&cu_TimerStop), "Failed to create stop event!"); 

    CUDA_CALL(cudaEventRecord(cu_TimerStart), "Failed to record start event!"); 
} 

float d_CUDATimerStop(void) 
{ 
    CUDA_CALL(cudaEventRecord(cu_TimerStop), "Failed to record stop event!"); 

    CUDA_CALL(cudaEventSynchronize(cu_TimerStop), "Failed to synch stop event!"); 

    float ms; 

    CUDA_CALL(cudaEventElapsedTime(&ms, cu_TimerStart, cu_TimerStop), "Failed to elapse events!"); 

    CUDA_CALL(cudaEventDestroy(cu_TimerStart), "Failed to destroy start event!"); 
    CUDA_CALL(cudaEventDestroy(cu_TimerStop), "Failed to destroy stop event!"); 

    return ms; 
} 

float* d_GetInv(float* L, int n) 
{ 
    cublasHandle_t cu_cublasHandle; 
    CUBLAS_CALL(cublasCreate(&cu_cublasHandle), "Failed to initialize cuBLAS!"); 

    float** adL; 
    float** adC; 
    float* dL; 
    float* dC; 
    int* dLUPivots; 
    int* dLUInfo; 

    size_t szA = n * n * sizeof(float); 

    CUDA_CALL(cudaMalloc(&adL, sizeof(float*)), "Failed to allocate adL!"); 
    CUDA_CALL(cudaMalloc(&adC, sizeof(float*)), "Failed to allocate adC!"); 
    CUDA_CALL(cudaMalloc(&dL, szA), "Failed to allocate dL!"); 
    CUDA_CALL(cudaMalloc(&dC, szA), "Failed to allocate dC!"); 
    CUDA_CALL(cudaMalloc(&dLUPivots, n * sizeof(int)), "Failed to allocate dLUPivots!"); 
    CUDA_CALL(cudaMalloc(&dLUInfo, sizeof(int)), "Failed to allocate dLUInfo!"); 

    CUDA_CALL(cudaMemcpy(dL, L, szA, cudaMemcpyHostToDevice), "Failed to copy to dL!"); 
    CUDA_CALL(cudaMemcpy(adL, &dL, sizeof(float*), cudaMemcpyHostToDevice), "Failed to copy to adL!"); 
    CUDA_CALL(cudaMemcpy(adC, &dC, sizeof(float*), cudaMemcpyHostToDevice), "Failed to copy to adC!"); 

    d_CUDATimerStart(); 

    CUBLAS_CALL(cublasSgetrfBatched(cu_cublasHandle, n, adL, n, dLUPivots, dLUInfo, 1), "Failed to perform LU decomp operation!"); 
    CUDA_CALL(cudaDeviceSynchronize(), "Failed to synchronize after kernel call!"); 

    CUBLAS_CALL(cublasSgetriBatched(cu_cublasHandle, n, (const float **)adL, n, dLUPivots, adC, n, dLUInfo, 1), "Failed to perform Inverse operation!"); 
    CUDA_CALL(cudaDeviceSynchronize(), "Failed to synchronize after kernel call!"); 

    float timed = d_CUDATimerStop(); 

    printf("cublas inverse in: %.5f ms.\n", timed); 

    float* res = (float*)malloc(szA); 

    CUDA_CALL(cudaMemcpy(res, dC, szA, cudaMemcpyDeviceToHost), "Failed to copy to res!"); 

    CUDA_CALL(cudaFree(adL), "Failed to free adL!"); 
    CUDA_CALL(cudaFree(adC), "Failed to free adC!"); 
    CUDA_CALL(cudaFree(dL), "Failed to free dL!"); 
    CUDA_CALL(cudaFree(dC), "Failed to free dC!"); 
    CUDA_CALL(cudaFree(dLUPivots), "Failed to free dLUPivots!"); 
    CUDA_CALL(cudaFree(dLUInfo), "Failed to free dLUInfo!"); 

    CUBLAS_CALL(cublasDestroy(cu_cublasHandle), "Failed to destroy cuBLAS!"); 

    return res; 
} 

int main() 
{ 
    int n = 1000; 
    float* L = (float*)malloc(n * n * sizeof(float)); 
    for(int i = 0; i < n * n; i++) 
     L[i] = ((float)rand()/(float)(RAND_MAX)); 

    float* inv = d_GetInv(L, n); 

    printf("done."); 
    _getch(); 

    return 0; 
} 

MATLAB code

A = rand(1000); 
tic 
X = inv(A); 
toc 

Informations système:

GPU: GTX 780 3Go

CPU: i7-4790S @ 3.20 GHz

Answer 1

Comme l'a dit @RobertCrovella, vous ne devriez pas utiliser les API de petite matrice pour une seule par lots grande inversion de matrice. Fondamentalement, vous pouvez utiliser la même méthode que dans votre code, mais avec la version non-lot de getrf() et getri() pour des performances maximales pour les grandes matrices. Pour getrf() vous pouvez le trouver ici.

http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-getrf

Pour getri(), bien que boîte à outils CUDA ne fournit pas getri() pour résoudre AX=I, où A est LU-facotored par getrf(), il fournit une getrs() pour résoudre AX=B. Tout ce que vous devez faire est de régler B=I avant d'appeler getrs().

http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-getrs