Как настроить аргументы функции cublas {t} symm()

Question

Эта функция выполняет симметричное умножение матричной матрицы с использованием CUDA. Хотя мне удалось использовать несимметричную версию «cublas {t} gemm()» Я не мог правильно использовать функцию «cublas {t} symm()».

Я знаю, что библиотека CUBLAS использует хранилище матричных столбцов. Я использую строчную матрицу C/C++, и я знаю, как решить эту проблему для «cublas {t} gemm()» путем замены входных матриц и т. Д. Однако я не мог решить ее для симметричного случая. Проблема в том, что даже если я использую хранилище матричных столбцов, я нахожу неожиданные результаты. Матрицы содержат комплексные поплавки (cuComplex). Я предполагаю, что у меня есть строковые матрицы. Вот код и выход:

// Matrix multiplication: C = A * B. 
// Host code. 
// 

// Utilities and system includes 
#include <assert.h> 
#include <helper_string.h> // helper for shared functions common to CUDA SDK samples 

// CUDA runtime 
#include <cuda_runtime.h> 
#include <cublas_v2.h> 

#ifndef min 
#define min(a,b) ((a < b) ? a : b) 
#endif 
#ifndef max 
#define max(a,b) ((a > b) ? a : b) 
#endif 

//////////////////////////////////////////////////////////////////////////////// 
// These are CUDA Helper functions (in addition to helper_cuda.h) 

void inline checkError(cublasStatus_t status, const char *msg) 
{ 
    if (status != CUBLAS_STATUS_SUCCESS) 
    { 
     printf("%s", msg); 
     exit(EXIT_FAILURE); 
    } 
} 
// end of CUDA Helper Functions 

// Allocates a matrix with random float entries. 
void randomCmplxInit(cuComplex *data, int size) 
{ 
    for (int i = 0; i < size; ++i) 
     data[i] = make_cuComplex(rand()/(float)RAND_MAX, rand()/(float)RAND_MAX); 
} 


//void initializeCUDA(int argc, char **argv, int &devID, int &iSizeMultiple, sMatrixSize &matrix_size) 
void initializeCUDA(int argc, char **argv, int &devID) 
{ 
    // By default, we use device 0, otherwise we override the device ID based on what is provided at the command line 
    cudaError_t error; 
    devID = 0; 
    int m,n,k; 

    if (checkCmdLineFlag(argc, (const char **)argv, "device")) 
    { 
     devID = getCmdLineArgumentInt(argc, (const char **)argv, "device"); 
     error = cudaSetDevice(devID); 

     if (error != cudaSuccess) 
     { 
      printf("cudaSetDevice returned error code %d, line(%d)\n", error, __LINE__); 
      exit(EXIT_FAILURE); 
     } 
    } 

    // get number of SMs on this GPU 
    error = cudaGetDevice(&devID); 
    cudaDeviceProp deviceProp; 
    error = cudaGetDeviceProperties(&deviceProp, devID); 

    printf("GPU Device %d: \"%s\" with compute capability %d.%d\n\n", devID, deviceProp.name, deviceProp.major, deviceProp.minor); 

    // use a larger block size for Fermi and above 
    int block_size = (deviceProp.major < 2) ? 16 : 32; 
} 

//////////////////////////////////////////////////////////////////////////////// 
//! Run a simple test matrix multiply using CUBLAS 
//////////////////////////////////////////////////////////////////////////////// 
int matrixMultiply(int argc, char **argv, int devID) 
{ 
    int i,j; 
    unsigned int m,n,k; 
    cudaDeviceProp deviceProp; 
    cudaError_t error; 

    error = cudaGetDeviceProperties(&deviceProp, devID); 

    if (error != cudaSuccess) 
    { 
     printf("cudaGetDeviceProperties returned error code %d, line(%d)\n", error, __LINE__); 
     exit(EXIT_FAILURE); 
    } 

    // use a larger block size for Fermi and above 
    int block_size = (deviceProp.major < 2) ? 16 : 32; 

    m=3; //number of rows of matrix op(A) and C. A--> (m x k) 
    n=2; //number of columns of matrix op(B) and C. B--> (k x n) 
    k=m; //number of columns of op(A) and rows of op(B). C--> (m x n) 

    // I want to compute C = A*B in row-major format, 
    //so I must find C(T)=B(T)A(T) = C(T)A in column-major format 

    // allocate host memory for matrices A and B 
    unsigned int size_A = m*(m+1)/2; //size of a symmetric matrix 
    unsigned int mem_size_A = sizeof(cuComplex) * size_A; 
    cuComplex *h_A = (cuComplex *)malloc(mem_size_A); 

    unsigned int size_B = m*n; 
    unsigned int mem_size_B = sizeof(cuComplex) * size_B; 
    cuComplex *h_B = (cuComplex *)malloc(mem_size_B); 

    // initialize host memory 
    for (i = 0; i < size_A; ++i) 
    h_A[i] = make_cuComplex((float)(i+1),(float)0); 

    for (i = 0; i < size_B; ++i) 
    h_B[i] = make_cuComplex((float)(i+2), (float)0); 

    // allocate device memory 
    cuComplex *d_A, *d_B, *d_C; 
    unsigned int size_C = m*n; 
    unsigned int mem_size_C = sizeof(cuComplex) * size_C; 

    // allocate host memory for the result 
    cuComplex *h_C  = (cuComplex *) malloc(mem_size_C); 
    cuComplex *h_CUBLAS = (cuComplex *) malloc(mem_size_C); 

    error = cudaMalloc((void **) &d_A, mem_size_A); 
    error = cudaMalloc((void **) &d_B, mem_size_B); 

    // copy host memory to device 
    error = cudaMemcpy(d_A, h_A, mem_size_A, cudaMemcpyHostToDevice); 
    error = cudaMemcpy(d_B, h_B, mem_size_B, cudaMemcpyHostToDevice); 
    error = cudaMalloc((void **) &d_C, mem_size_C); 

    // setup execution parameters 
    dim3 threads(block_size, block_size); 
    dim3 grid(n/threads.x, m/threads.y); 

    // create and start timer 
    printf("Computing result using CUBLAS..."); 

    // CUBLAS version 2.0 
    { 
     cublasHandle_t handle; 
     cublasStatus_t ret; 

     ret = cublasCreate(&handle); 

     if (ret != CUBLAS_STATUS_SUCCESS) 
     { 
      printf("cublasCreate returned error code %d, line(%d)\n", ret, __LINE__); 
      exit(EXIT_FAILURE); 
     } 

     const cuComplex alpha = make_cuComplex(1.0f,0.0f); 
     const cuComplex beta = make_cuComplex(0.0f,0.0f); 
     //Perform operation with cublas 
     ret = cublasCsymm(handle, CUBLAS_SIDE_RIGHT, CUBLAS_FILL_MODE_UPPER, n,m,&alpha,d_A,m,d_B,m,&beta,d_C,m); 

     // copy result from device to host 
     error = cudaMemcpy(h_CUBLAS, d_C, mem_size_C, cudaMemcpyDeviceToHost); 

     checkError(cublasDestroy(handle), "cublasDestroy() error!\n"); 
    } 


    printf ("\nComputations completed.\n\n"); 

    printf (" symm matrix A: \n"); 
    int s=0; 
    for (i=0; i<min(m,4); i++) { 
     for (j=0; j<=i; j++) { 
     //printf ("%7.5G + j(%7.5G)", h_A[j+i*k].x,h_A[j+i*k].y); 
     printf ("%7.5G", h_A[s].x); 
     s++; 
     } 
     printf ("\n"); 
    } 

    printf ("\n matrix B: \n"); 
    for (i=0; i<min(k,4); i++) { 
     for (j=0; j<min(n,4); j++) { 
     //printf ("%7.5G + j(%7.5G)", h_B[j+i*n].x,h_B[j+i*n].y); 
      printf ("%7.5G", h_B[j+i*n].x); 
     } 
     printf ("\n"); 
    } 

    printf ("\n matrix C=A*B: \n"); 
    for (i=0; i<min(m,4); i++) { 
     for (j=0; j<min(n,4); j++) { 
     //printf ("%7.5G + j(%7.5G)", h_CUBLAS[j+i*n].x,h_CUBLAS[j+i*n].y); 
      printf ("%7.5G", h_CUBLAS[j+i*n].x); 
     } 
     printf ("\n"); 
    } 

    // clean up memory 
    free(h_A); 
    free(h_B); 
    free(h_C); 
    //free(reference); 
    cudaFree(d_A); 
    cudaFree(d_B); 
    cudaFree(d_C); 

    cudaDeviceReset(); 
} 

//////////////////////////////////////////////////////////////////////////////// 
// Program main 
//////////////////////////////////////////////////////////////////////////////// 
int main(int argc, char **argv) 
{ 
    printf("[Matrix Multiply CUBLAS] - Starting...\n"); 
    int devID = 0, sizeMult = 5; 
    initializeCUDA(argc, argv, devID); 
    int matrix_result = matrixMultiply(argc, argv, devID); 
} 

Я полагаю, что у меня есть следующие матрицы для умножения:

A =

1  2  4 
2  3  5 
4  5  6 

B =

2  3 
4  5 
6  7 

и рассчитывать на получение

A * B =

34 41 
46 56 
64 79 

Но полученный ВЫВОД выглядит следующим образом:

symm matrix A:  
    1  
    2  3 
    4  5  6 

matrix B: 
    2  3 
    4  5 
    6  7 

matrix C=A*B: 
    78 90 
    74 97 
    114 146 

Что я упускаю в этом коде? Вероятно, аргументы функции «cublasCsymm» неверны.

Спасибо, Каган

Answer 1

EDIT: Основываясь на поставленные вопросы ниже, я решил вновь работать мой ответ и пример кода. Вы можете обрабатывать хранилище больших строк без транспонирования, по крайней мере, для этих операций. И это наблюдение дополнительно облегчается тем фактом, что функция symm не использует упакованное хранилище.

Так, чтобы ответить на дополнительные вопросы:

1. функция cublasCsymm не использует упакованный формат хранения (например, некоторые другие функции, такие как cublasCspmv, например), так как функция cublasCsymm предназначена для дублирования функциональности соответствующий netlib function, который также не использует формат упакованного хранилища. Основываясь на моем обзоре API cublas, я не вижу функцию умножения матричной матрицы с симметричной упаковкой.
2. Может использовать с использованием хранилища строк (например, C-style) с помощью cublas, без переноса, по крайней мере для этих операций (матрица-матрица размножается, без упакованного хранения), следуя рекомендациям here.

Ниже приведена переработанная версия моего предыдущего примера, которая включает информацию в пункте 2 выше.

// Matrix multiplication: C = A * B. 
// Host code. 
// 

// Utilities and system includes 
#include <assert.h> 
#include <helper_string.h> // helper for shared functions common to CUDA SDK sa 
mples 

// CUDA runtime 
#include <cuda_runtime.h> 
#include <cublas_v2.h> 

// error check macros 
#define cudaCheckErrors(msg) \ 
    do { \ 
     cudaError_t __err = cudaGetLastError(); \ 
     if (__err != cudaSuccess) { \ 
      fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \ 
       msg, cudaGetErrorString(__err), \ 
       __FILE__, __LINE__); \ 
      fprintf(stderr, "*** FAILED - ABORTING\n"); \ 
      exit(1); \ 
     } \ 
    } while (0) 

// for CUBLAS V2 API 
#define cublasCheckErrors(fn) \ 
    do { \ 
     cublasStatus_t __err = fn; \ 
     if (__err != CUBLAS_STATUS_SUCCESS) { \ 
      fprintf(stderr, "Fatal cublas error: %d (at %s:%d)\n", \ 
       (int)(__err), \ 
       __FILE__, __LINE__); \ 
      fprintf(stderr, "*** FAILED - ABORTING\n"); \ 
      exit(1); \ 
     } \ 
    } while (0) 



#ifndef min 
#define min(a,b) ((a < b) ? a : b) 
#endif 
#ifndef max 
#define max(a,b) ((a > b) ? a : b) 
#endif 

//////////////////////////////////////////////////////////////////////////////// 
// These are CUDA Helper functions (in addition to helper_cuda.h) 

void inline checkError(cublasStatus_t status, const char *msg) 
{ 
    if (status != CUBLAS_STATUS_SUCCESS) 
    { 
     printf("%s", msg); 
     exit(EXIT_FAILURE); 
    } 
} 
// end of CUDA Helper Functions 

// Allocates a matrix with random float entries. 
void randomCmplxInit(cuComplex *data, int size) 
{ 
    for (int i = 0; i < size; ++i) 
     data[i] = make_cuComplex(rand()/(float)RAND_MAX, rand()/(float)RAND 
_MAX); 
} 


//void initializeCUDA(int argc, char **argv, int &devID, int &iSizeMultiple, sMa 
trixSize &matrix_size) 
void initializeCUDA(int argc, char **argv, int &devID) 
{ 
    // By default, we use device 0, otherwise we override the device ID based on 
what is provided at the command line 
    cudaError_t error; 
    devID = 0; 

    if (checkCmdLineFlag(argc, (const char **)argv, "device")) 
    { 
     devID = getCmdLineArgumentInt(argc, (const char **)argv, "device"); 
     error = cudaSetDevice(devID); 

     if (error != cudaSuccess) 
     { 
      printf("cudaSetDevice returned error code %d, line(%d)\n", error, __ 
LINE__); 
      exit(EXIT_FAILURE); 
     } 
    } 

    // get number of SMs on this GPU 
    error = cudaGetDevice(&devID); 
    cudaDeviceProp deviceProp; 
    error = cudaGetDeviceProperties(&deviceProp, devID); 

    printf("GPU Device %d: \"%s\" with compute capability %d.%d\n\n", devID, dev 
iceProp.name, deviceProp.major, deviceProp.minor); 

} 


//////////////////////////////////////////////////////////////////////////////// 
//! Run a simple test matrix multiply using CUBLAS 
//////////////////////////////////////////////////////////////////////////////// 
int matrixMultiply(int argc, char **argv, int devID) 
{ 
    int i,j; 
    unsigned int m,n,k; 
    cudaDeviceProp deviceProp; 
    cudaError_t error; 

    error = cudaGetDeviceProperties(&deviceProp, devID); 

    if (error != cudaSuccess) 
    { 
     printf("cudaGetDeviceProperties returned error code %d, line(%d)\n", error, __LINE__); 
     exit(EXIT_FAILURE); 
    } 

    // use a larger block size for Fermi and above 

    m=3; //number of rows of matrix op(A) and C. A--> (m x k) 
    n=2; //number of columns of matrix op(B) and C. B--> (k x n) 
    k=m; //number of columns of op(A) and rows of op(B). C--> (m x n) 

    // I want to compute C = A*B in row-major format, 
    //so I must find C(T)=B(T)A(T) = C(T)A in column-major format 

    // allocate host memory for matrices A and B 
    unsigned int size_A = m*m; //size of a symmetric matrix 
    printf("size_A = %d\n", size_A); 
    unsigned int mem_size_A = sizeof(cuComplex) * size_A; 
    cuComplex *h_A = (cuComplex *)malloc(mem_size_A); 

    unsigned int size_B = m*n; 
    unsigned int mem_size_B = sizeof(cuComplex) * size_B; 
    cuComplex *h_B = (cuComplex *)malloc(mem_size_B); 

    // initialize host memory 
// for (i = 0; i < size_A; ++i) 
// h_A[i] = make_cuComplex((float)(i+1),(float)0); 
    h_A[0] = make_cuComplex((float)1, (float)0); 
    h_A[1] = make_cuComplex((float)2, (float)0); 
    h_A[2] = make_cuComplex((float)4, (float)0); 
    h_A[3] = make_cuComplex((float)0, (float)0); 
    h_A[4] = make_cuComplex((float)3, (float)0); 
    h_A[5] = make_cuComplex((float)5, (float)0); 
    h_A[6] = make_cuComplex((float)0, (float)0); 
    h_A[7] = make_cuComplex((float)0, (float)0); 
    h_A[8] = make_cuComplex((float)6, (float)0); 

// for (i = 0; i < size_B; ++i) 
// h_B[i] = make_cuComplex((float)(i+2), (float)0); 
    h_B[0] = make_cuComplex((float)2, (float)0); 
    h_B[1] = make_cuComplex((float)3, (float)0); 
    h_B[2] = make_cuComplex((float)4, (float)0); 
    h_B[3] = make_cuComplex((float)5, (float)0); 
    h_B[4] = make_cuComplex((float)6, (float)0); 
    h_B[5] = make_cuComplex((float)7, (float)0); 

    // allocate device memory 
    cuComplex *d_A, *d_B, *d_C; 
    unsigned int size_C = m*n; 
    unsigned int mem_size_C = sizeof(cuComplex) * size_C; 

    // allocate host memory for the result 
    cuComplex *h_C  = (cuComplex *) malloc(mem_size_C); 
    cuComplex *h_CUBLAS = (cuComplex *) malloc(mem_size_C); 

    error = cudaMalloc((void **) &d_A, mem_size_A); 
    error = cudaMalloc((void **) &d_B, mem_size_B); 

    // copy host memory to device 
    error = cudaMemcpy(d_A, h_A, mem_size_A, cudaMemcpyHostToDevice); 
    error = cudaMemcpy(d_B, h_B, mem_size_B, cudaMemcpyHostToDevice); 
    error = cudaMalloc((void **) &d_C, mem_size_C); 


    // create and start timer 
    printf("Computing result using CUBLAS..."); 

    // CUBLAS version 2.0 
    { 
     cublasHandle_t handle; 
     cublasStatus_t ret; 

     ret = cublasCreate(&handle); 

     if (ret != CUBLAS_STATUS_SUCCESS) 
     { 
      printf("cublasCreate returned error code %d, line(%d)\n", ret, __LINE__); 
      exit(EXIT_FAILURE); 
     } 

     const cuComplex alpha = make_cuComplex(1.0f,0.0f); 
     const cuComplex beta = make_cuComplex(0.0f,0.0f); 
     //Perform operation with cublas 
     ret = cublasCsymm(handle, CUBLAS_SIDE_RIGHT, CUBLAS_FILL_MODE_LOWER, n,m,&alpha,d_A,m,d_B,n,&beta,d_C,n); 

     if (ret != CUBLAS_STATUS_SUCCESS) 
     { 
      printf("cublasCsymm returned error code %d, line(%d)\n", ret, __LINE__); 
      exit(EXIT_FAILURE); 
     } 

     // copy result from device to host 
     error = cudaMemcpy(h_CUBLAS, d_C, mem_size_C, cudaMemcpyDeviceToHost); 

     checkError(cublasDestroy(handle), "cublasDestroy() error!\n"); 
    } 


    printf ("\nComputations completed.\n\n"); 

    printf (" symm matrix A: \n"); 
// int s=0; 
    for (i=0; i<min(m,4); i++) { 
     for (j=0; j<min(m,4); j++) { 
     //printf ("%7.5G + j(%7.5G)", h_A[j+i*k].x,h_A[j+i*k].y); 
//  printf ("%7.5G", h_A[s].x); 
     printf ("%7.5G", h_A[j+(i*m)].x); 
//  s++; 
     } 
     printf ("\n"); 
    } 

    printf ("\n matrix B: \n"); 
    for (i=0; i<min(k,4); i++) { 
     for (j=0; j<min(n,4); j++) { 
     //printf ("%7.5G + j(%7.5G)", h_B[j+i*n].x,h_B[j+i*n].y); 
      printf ("%7.5G", h_B[j+(i*n)].x); 
     } 
     printf ("\n"); 
    } 

    printf ("\n matrix C=A*B: \n"); 
    for (i=0; i<min(m,4); i++) { 
     for (j=0; j<min(n,4); j++) { 
     //printf ("%7.5G + j(%7.5G)", h_CUBLAS[j+i*n].x,h_CUBLAS[j+i*n].y); 
      printf ("%7.5G", h_CUBLAS[j+(i*n)].x); 
     } 
     printf ("\n"); 
    } 

    // clean up memory 
    free(h_A); 
    free(h_B); 
    free(h_C); 
    //free(reference); 
    cudaFree(d_A); 
    cudaFree(d_B); 
    cudaFree(d_C); 

    cudaDeviceReset(); 
    return 0; 
} 

//////////////////////////////////////////////////////////////////////////////// 
// Program main 
//////////////////////////////////////////////////////////////////////////////// 
int main(int argc, char **argv) 
{ 
    printf("[Matrix Multiply CUBLAS] - Starting...\n"); 
    int devID = 0; 
    initializeCUDA(argc, argv, devID); 
    int matrix_result = matrixMultiply(argc, argv, devID); 
    cudaCheckErrors("some error"); 
    return 0; 
} 
$ ./t213 
[Matrix Multiply CUBLAS] - Starting... 
GPU Device 0: "Tesla M2070" with compute capability 2.0 

size_A = 9 
Computing result using CUBLAS... 
Computations completed. 

symm matrix A: 
     1  2  4 
     0  3  5 
     0  0  6 

matrix B: 
     2  3 
     4  5 
     6  7 

matrix C=A*B: 
    34  41 
    46  56 
    64  79 
$ 

ORIGINAL РЕПЛИКА:

Несколько проблем:

• Когда я запускаю свой код, как вы его размещены прямо сейчас, я не получаю результаты, которые вы показываете. Вот что я получаю:

  [Matrix Multiply CUBLAS] - Starting... 
  GPU Device 0: "Tesla M2070" with compute capability 2.0 
  
  Computing result using CUBLAS... 
  Computations completed. 
  
  symm matrix A: 
      1 
      2  3 
      4  5  6 
  
  matrix B: 
      2  3 
      4  5 
      6  7 
  
  matrix C=A*B: 
      -131 -128 
      260 -122 
      -115 266 
  

Код компилируется с количеством предупреждений, а также вы не делаете надлежащую проверку ошибок (например, вы не проверяя возвращаемое значение cublasCsymm

• вы желаете, чтобы умножить C = A * B Это означает, что а на СЛЕВА, но вы передаете CUBLAS_SIDE_RIGHT в cublasCsymm Несколько других cublasCsymm параметры были неправильными. Я думаю, может быть, вы думали, что можете сделать A*B как (B (T) * A (T)), но это работает только для квадратных матриц. Не уверен, что вы думаете, точно.

• У вас есть хранилище строк на ваших матрицах и передача их в кубы, которые интерпретируют их в основном порядке. Для следующей матрицы:

  1 2 
  3 4 
  

строка-мажорного хранение выглядит следующим образом:

1 2 3 4 

столбцы хранение выглядит следующим образом:

1 3 2 4 

Вы можете перенести эти матрицы если хотите, используя cublasCgeam или вы можете вручную изменить свое хранилище.

• Вы делаете какое-то предположение о каком-то формат сжатого для хранения симметричной матрицы A, которая не является правильной. Прочтите внимательно определение storage type. Он не говорит, что часть матрицы, «поставленная» или «присутствует», говорит, что часть матрицы, которая находится , заполнена.

Вот полный код, который имеет вышеуказанные проблемы исправлено:

// Matrix multiplication: C = A * B. 
// Host code. 
// 

// Utilities and system includes 
#include <assert.h> 
#include <helper_string.h> // helper for shared functions common to CUDA SDK sa 
mples 

// CUDA runtime 
#include <cuda_runtime.h> 
#include <cublas_v2.h> 

// error check macros 
#define cudaCheckErrors(msg) \ 
    do { \ 
     cudaError_t __err = cudaGetLastError(); \ 
     if (__err != cudaSuccess) { \ 
      fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \ 
       msg, cudaGetErrorString(__err), \ 
       __FILE__, __LINE__); \ 
      fprintf(stderr, "*** FAILED - ABORTING\n"); \ 
      exit(1); \ 
     } \ 
    } while (0) 

// for CUBLAS V2 API 
#define cublasCheckErrors(fn) \ 
    do { \ 
     cublasStatus_t __err = fn; \ 
     if (__err != CUBLAS_STATUS_SUCCESS) { \ 
      fprintf(stderr, "Fatal cublas error: %d (at %s:%d)\n", \ 
       (int)(__err), \ 
       __FILE__, __LINE__); \ 
      fprintf(stderr, "*** FAILED - ABORTING\n"); \ 
      exit(1); \ 
     } \ 
    } while (0) 



#ifndef min 
#define min(a,b) ((a < b) ? a : b) 
#endif 
#ifndef max 
#define max(a,b) ((a > b) ? a : b) 
#endif 

//////////////////////////////////////////////////////////////////////////////// 
// These are CUDA Helper functions (in addition to helper_cuda.h) 


void inline checkError(cublasStatus_t status, const char *msg) 
{ 
    if (status != CUBLAS_STATUS_SUCCESS) 
    { 
     printf("%s", msg); 
     exit(EXIT_FAILURE); 
    } 
} 
// end of CUDA Helper Functions 

// Allocates a matrix with random float entries. 
void randomCmplxInit(cuComplex *data, int size) 
{ 
    for (int i = 0; i < size; ++i) 
     data[i] = make_cuComplex(rand()/(float)RAND_MAX, rand()/(float)RAND_MAX); 
} 


//void initializeCUDA(int argc, char **argv, int &devID, int &iSizeMultiple, sMatrixSize &matrix_size) 
void initializeCUDA(int argc, char **argv, int &devID) 
{ 
    // By default, we use device 0, otherwise we override the device ID based on what is provided at the command line 
    cudaError_t error; 
    devID = 0; 

    if (checkCmdLineFlag(argc, (const char **)argv, "device")) 
    { 
     devID = getCmdLineArgumentInt(argc, (const char **)argv, "device"); 
     error = cudaSetDevice(devID); 

     if (error != cudaSuccess) 
     { 
      printf("cudaSetDevice returned error code %d, line(%d)\n", error, __LINE__); 
      exit(EXIT_FAILURE); 
     } 
    } 

    // get number of SMs on this GPU 
    error = cudaGetDevice(&devID); 
    cudaDeviceProp deviceProp; 
    error = cudaGetDeviceProperties(&deviceProp, devID); 

    printf("GPU Device %d: \"%s\" with compute capability %d.%d\n\n", devID, deviceProp.name, deviceProp.major, deviceProp.minor); 

} 

//////////////////////////////////////////////////////////////////////////////// 
//! Run a simple test matrix multiply using CUBLAS 
//////////////////////////////////////////////////////////////////////////////// 
int matrixMultiply(int argc, char **argv, int devID) 
{ 
    int i,j; 
    unsigned int m,n,k; 
    cudaDeviceProp deviceProp; 
    cudaError_t error; 

    error = cudaGetDeviceProperties(&deviceProp, devID); 

    if (error != cudaSuccess) 
    { 
     printf("cudaGetDeviceProperties returned error code %d, line(%d)\n", error, __LINE__); 
     exit(EXIT_FAILURE); 
    } 

    // use a larger block size for Fermi and above 

    m=3; //number of rows of matrix op(A) and C. A--> (m x k) 
    n=2; //number of columns of matrix op(B) and C. B--> (k x n) 
    k=m; //number of columns of op(A) and rows of op(B). C--> (m x n) 

    // I want to compute C = A*B in row-major format, 
    //so I must find C(T)=B(T)A(T) = C(T)A in column-major format 

    // allocate host memory for matrices A and B 
    unsigned int size_A = m*m; //size of a symmetric matrix 
    printf("size_A = %d\n", size_A); 
    unsigned int mem_size_A = sizeof(cuComplex) * size_A; 
    cuComplex *h_A = (cuComplex *)malloc(mem_size_A); 

    unsigned int size_B = m*n; 
    unsigned int mem_size_B = sizeof(cuComplex) * size_B; 
    cuComplex *h_B = (cuComplex *)malloc(mem_size_B); 

    // initialize host memory 
// for (i = 0; i < size_A; ++i) 
// h_A[i] = make_cuComplex((float)(i+1),(float)0); 
    h_A[0] = make_cuComplex((float)1, (float)0); 
    h_A[1] = make_cuComplex((float)2, (float)0); 
    h_A[2] = make_cuComplex((float)4, (float)0); 
    h_A[3] = make_cuComplex((float)0, (float)0); 
    h_A[4] = make_cuComplex((float)3, (float)0); 
    h_A[5] = make_cuComplex((float)5, (float)0); 
    h_A[6] = make_cuComplex((float)0, (float)0); 
    h_A[7] = make_cuComplex((float)0, (float)0); 
    h_A[8] = make_cuComplex((float)6, (float)0); 

// for (i = 0; i < size_B; ++i) 
// h_B[i] = make_cuComplex((float)(i+2), (float)0); 
    h_B[0] = make_cuComplex((float)2, (float)0); 
    h_B[1] = make_cuComplex((float)4, (float)0); 
    h_B[2] = make_cuComplex((float)6, (float)0); 
    h_B[3] = make_cuComplex((float)3, (float)0); 
    h_B[4] = make_cuComplex((float)5, (float)0); 
    h_B[5] = make_cuComplex((float)7, (float)0); 

    // allocate device memory 
    cuComplex *d_A, *d_B, *d_C; 
    unsigned int size_C = m*n; 
    unsigned int mem_size_C = sizeof(cuComplex) * size_C; 

    // allocate host memory for the result 
    cuComplex *h_C  = (cuComplex *) malloc(mem_size_C); 
    cuComplex *h_CUBLAS = (cuComplex *) malloc(mem_size_C); 

    error = cudaMalloc((void **) &d_A, mem_size_A); 
    error = cudaMalloc((void **) &d_B, mem_size_B); 

    // copy host memory to device 
    error = cudaMemcpy(d_A, h_A, mem_size_A, cudaMemcpyHostToDevice); 
    error = cudaMemcpy(d_B, h_B, mem_size_B, cudaMemcpyHostToDevice); 
    error = cudaMalloc((void **) &d_C, mem_size_C); 


    // create and start timer 
    printf("Computing result using CUBLAS..."); 

    // CUBLAS version 2.0 
    { 
     cublasHandle_t handle; 
     cublasStatus_t ret; 

     ret = cublasCreate(&handle); 

     if (ret != CUBLAS_STATUS_SUCCESS) 
     { 
      printf("cublasCreate returned error code %d, line(%d)\n", ret, __LINE__); 
      exit(EXIT_FAILURE); 
     } 

     const cuComplex alpha = make_cuComplex(1.0f,0.0f); 
     const cuComplex beta = make_cuComplex(0.0f,0.0f); 
     //Perform operation with cublas 
     ret = cublasCsymm(handle, CUBLAS_SIDE_LEFT, CUBLAS_FILL_MODE_LOWER, m,n,&alpha,d_A,m,d_B,m,&beta,d_C,m); 
     if (ret != CUBLAS_STATUS_SUCCESS) 
     { 
      printf("cublasCsymm returned error code %d, line(%d)\n", ret, __LINE__); 
      exit(EXIT_FAILURE); 
     } 

Вот результат:

[Matrix Multiply CUBLAS] - Starting... 
GPU Device 0: "Tesla M2070" with compute capability 2.0 

size_A = 9 
Computing result using CUBLAS... 
Computations completed. 

symm matrix A: 
     1  0  0 
     2  3  0 
     4  5  6 

matrix B: 
     2  3 
     4  5 
     6  7 

matrix C=A*B: 
    34  41 
    46  56 
    64  79