cublasDGemm résultats étranges

Question

J'essaie d'utiliser cublasDgemm() pour calculer le produit d'une matrice par sa transposition. La matrice d'entrée et la sortie j'attends de mon code sont les suivants (A et C respectivement):

| 1 4 7 |  | 66 78 | 
A = | 2 5 8 | C = | 78 93 | 

Cependant, j'obtiens des résultats étranges et il est un peu difficile pour moi de comprendre les dimensions cublas/cuda utilise (colonne majeure). Tout conseil serait apprécié!

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <cuda_runtime.h> 
#include "cublas_v2.h" 
#define M 3 
#define N 2 
#define IDX2C(i,j,ld) (((j)*(ld))+(i)) 

int main (void){ 
    cudaError_t cudaStat;  
    cublasStatus_t stat; 
    cublasHandle_t handle; 
    int i, j; 
    double *devPtrA, *devPtrC; 
    double *a = 0, *c = 0; 

    const double alpha = 1; 
    const double beta = 0; 

    // initialize host arrays 
    a = (double *)malloc (M * N * sizeof (*a)); 
    c = (double *)malloc (N * N * sizeof (*c)); 
    if (!a || !c) { 
     printf ("host memory allocation failed"); 
     return EXIT_FAILURE; 
    } 

    // fill input array 
    for (j = 0; j < N; j++) { 
     for (i = 0; i < M; i++) { 
      a[IDX2C(i,j,M)] = (double)(i * M + j + 1); 
      printf ("%7.0f", a[IDX2C(i,j,M)]); 
     } 
     printf ("\n"); 
    } 

    // set device to 0 (for double processing) 
    cudaStat = cudaSetDevice(0); 
    if (cudaStat != cudaSuccess) { 
     printf("could not set device 0"); 
     return EXIT_FAILURE; 
    } 

    // allocate device arrays 
    cudaStat = cudaMalloc ((void**)&devPtrA, M*N*sizeof(*a)); 
    if (cudaStat != cudaSuccess) { 
     printf ("device memory allocation of A failed"); 
     return EXIT_FAILURE; 
    } 
    cudaStat = cudaMalloc ((void**)&devPtrC, N*N*sizeof(*c)); 
    if (cudaStat != cudaSuccess) { 
     printf ("device memory allocation of C failed"); 
     return EXIT_FAILURE; 
    } 

    // create the cublas handle 
    stat = cublasCreate(&handle); 
    if (stat != CUBLAS_STATUS_SUCCESS) { 
     printf ("CUBLAS initialization failed\n"); 
     return EXIT_FAILURE; 
    } 

    // set the matrix a 
    stat = cublasSetMatrix (M, N, sizeof(*a), a, M, devPtrA, M); 
    if (stat != CUBLAS_STATUS_SUCCESS) { 
     printf ("data download failed"); 
     cudaFree (devPtrA); 
     cudaFree (devPtrC); 
     cublasDestroy(handle); 
     return EXIT_FAILURE; 
    } 

    stat = cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, M, M, N, &alpha, devPtrA, M, devPtrA, M, &beta, devPtrC, M); 
    if (stat!= CUBLAS_STATUS_SUCCESS) { 
     switch (stat) { 
      case CUBLAS_STATUS_NOT_INITIALIZED: 
       printf("CUBLAS_STATUS_NOT_INITIALIZED\n"); 
       break; 
      case CUBLAS_STATUS_INVALID_VALUE: 
       printf("CUBLAS_STATUS_INVALID_VALUE\n"); 
       break; 
      case CUBLAS_STATUS_ARCH_MISMATCH: 
       printf("CUBLAS_STATUS_ARCH_MISMATCH\n"); 
       break; 
      case CUBLAS_STATUS_EXECUTION_FAILED: 
       printf("CUBLAS_STATUS_EXECUTION_FAILED\n"); 
       break; 
      default: 
       printf("??\n"); 
     } 

     printf("Error: %d\n", (int)stat); 
     cudaFree (devPtrA); 
     cudaFree (devPtrC); 
     cublasDestroy(handle); 
     return EXIT_FAILURE; 
    } 

    // get matrix c 
    stat = cublasGetMatrix (N, N, sizeof(*c), devPtrC, N, c, N); 
    if (stat != CUBLAS_STATUS_SUCCESS) { 
     printf ("data upload failed"); 
     cudaFree (devPtrC); 
     cublasDestroy(handle); 
     return EXIT_FAILURE; 
    } 

    // cleanup cuda/cublas 
    cudaFree (devPtrA); 
    cudaFree (devPtrC); 
    cublasDestroy(handle); 

    // print result 
    for (j = 0; j < N; j++) { 
     for (i = 0; i < N; i++) { 
      printf ("%7.0f", c[IDX2C(i,j,M)]); 
     } 
     printf ("\n"); 
    } 

    // clear host data 
    free(a); 
    free(c); 
    return EXIT_SUCCESS; 
} 

Answer 1

Le premier problème est que vous remplissez votre matrice A au format rangée-majeur. Pour résoudre ce problème, il suffit d'échanger les indices i et j. Dans le format de colonne principale, la dimension de premier plan devrait être le nombre de lignes, dans votre cas, N.

for (j = 0; j < N; j++) { 
    for (i = 0; i < M; i++) { 
     a[IDX2C(j,i,N)] = (double)(i * M + j + 1); 
     printf ("%7.0f", a[IDX2C(j,i,N)]); 
    } 
    printf ("\n"); 
} 

Vous échangez aussi des dimensions dans l'appel cublasDgemm, il devrait ressembler à:

stat = cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, N, N, M, &alpha, devPtrA, N, devPtrA, N, &beta, devPtrC, N); 

Et à la fin, vous avez utilisé M comme la dimension principale de la matrice C, où il aurait dû être N:

printf ("%7.0f", c[IDX2C(i,j,N)]);