2013-08-29 2 views
1

Я ранее разместил вопрос о умножении матричных векторов в CUDA и о написании собственного ядра. Сделав это, я решил реализовать свою проблему с помощью CUBLAS, как это было предложено некоторыми пользователями (спасибо @Robert Crovella) на SO в надежде на достижение более высокой производительности (мой проект зависит от производительности).CUDA/CUBLAS Матрица-векторное умножение

Просто уточнить: я хочу умножить матрицу NxN с вектором 1xN.

Я смотрел на код, наклеенный ниже на пару дней, и я не могу понять, почему умножение дает мне неправильный результат. Я боюсь, что я вызываю проблемы, используя < vector> массивы (это часть гораздо более крупной системы, использующей эти типы данных). Я не хочу использовать этот поток как инструмент отладки, но я думаю, что это также будет полезно для других пользователей, пытающихся достичь этого, поскольку я не сталкивался с особенно сложным источником в Интернете по моей конкретной проблеме (и для куб. v2 API). Заранее спасибо!

#include <cuda.h> 
#include <vector> 
#include <iostream> 
#include <fstream> 
#include <stdio.h> 
#include <stdlib.h> 
#include <cmath> 
#include <cublas_v2.h> 
#include <time.h> 

//#include "timenow.cu" 

// 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) 

// random data filler 
void fillvector(float *data, int N){ 
    for(int i=0; i<N; i++){ 
     data[i] = float(rand() % 10); 
    } 
} 

//printer 
void printer(bool printOut, float *data, int N){ 
    if(printOut == true){ 
    for(int i=0; i<N; i++){ 
     printf("%2.1f ", data[i]); 
    } 
    printf("\n"); 
    } 
} 

///////////////////////////////////////////////////////////////////// 
///////////////////////////////////////////////////////////////////// 

int main(){ 

bool printOut = true; 

int N; 
std::cout << "Enter N: " ; 
std::cin >> N; 

std::vector<float> x0; 
x0.resize(N); 

std::vector<float> p; 
p.resize(N*N); 

// matrix A 
std::vector<float> A[N]; 
for(int i=0;i<N;i++){ 
     A[i].resize(N); 
    fillvector(A[i].data(), N); 
    printer(printOut, A[i].data(), N); 
} 
printf("\n"); 
fillvector(x0.data(), N); 
printer(printOut, x0.data(), N); 

printf("\nStarting CUDA computation..."); 
///double startTime = timenow(); 

// device pointers 
float *d_A, *d_p, *d_b, *d_x0, *d_v, *d_temp; 

cudaMalloc((void**)&d_A, N*N*sizeof(float)); 
cudaMalloc((void**)&d_temp, N*sizeof(float)); 
cudaMalloc((void**)&d_x0, N*sizeof(float)); 
cudaCheckErrors("cuda malloc fail"); 

// might need to flatten A... 
cublasSetVector(N, sizeof(float), &x0, 1, d_x0, 1); 
//daMemcpy(d_x0, &x0, N*sizeof(float), cudaMemcpyHostToDevice); 
cublasSetMatrix(N, N, sizeof(float), &A, N, d_A, N); 
cudaCheckErrors("cuda memcpy of A or x0 fail"); 

float *temp; 
temp = (float *)malloc(N*sizeof(temp)); 

cublasHandle_t handle; 
cublasCheckErrors(cublasCreate(&handle)); 

float alpha = 1.0f; 
float beta = 0.0f; 
cublasCheckErrors(cublasSgemv(handle, CUBLAS_OP_N, N, N, &alpha, d_A, N, d_x0, 1, &beta, d_temp, 1)); 

cublasGetVector(N, sizeof(float), &temp, 1, d_temp, 1); 
//cudaMemcpy(temp, d_temp, N*sizeof(float), cudaMemcpyDeviceToHost); 
cudaCheckErrors("returning to host failed"); 

printf("\n"); 
printer(printOut, temp, N); 

/*alpha = -1.0; 
cublasSaxpy(handle, N, &alpha, d_temp, 1, d_v, 1); 
cublasGetVector(N, sizeof(float) * N, d_v, 1, &v, 1); 
printf("\n"); 
for(int i=0; i<N; i++){ 
    printf("%2.1f ",v[i]); 
}*/ 

printf("\nFinished CUDA computations..."); 
//double endTime = timenow(); 

//double timeDiff = endTime - startTime; 
//printf("\nRuntime: %2.3f seconds \n", timeDiff); 

cudaFree(d_temp); 
cudaFree(d_A); 
cudaFree(d_p); 
cudaFree(d_x0); 

return 0; 
} 
+0

Что это ваш вопрос? Если это умножение дает неверный результат, пожалуйста, дайте нам представление о результатах, которые вы получаете, и о результатах, которые вы ожидаете. Кроме того, вы не выполняете проверку ошибок cublas при получении/наборе матричных/векторных вызовов. –

+0

Ваша программа, похоже, допускает любой ввод. Когда я ввожу N = 2, я получаю сообщение об ошибке cuda с недопустимым аргументом в вашем сообщении «возврат к хосту». В чем ваш вопрос? –

+0

Привет, Роберт, я хотел бы умножить квадратную матрицу NxN на вектор 1xN. N должен быть любого размера. – rex

ответ

2
  • Мы не ссылаться на первый элемент вектора таким образом:

    cublasSetVector(N, sizeof(float), &x0, 1, d_x0, 
    

Вместо этого вы должны сделать это:

cublasSetVector(N, sizeof(float), &(x0[0]), 1, d_x0, 1); 

и также для ваш SetMatrix ссылка на вызов A:

cublasSetMatrix(N, N, sizeof(float), &(A[0]), N, d_A, N); 
  • Ваш GetVector вызов имеет 2 ошибки:

    cublasGetVector(N, sizeof(float), &temp, 1, d_temp, 1); 
    

У вас есть temp и d_temp параметры reversed (вы копируете от устройства к хосту), и вы не должны брать адрес temp: это уже указатель. Так что это:

cublasGetVector(N, sizeof(float), d_temp, 1, temp, 1); 
  • Вы не делаете надлежащую проверку ошибок на всех cublas вызовов, такие, как ваши получить/установить матрицу/вектор вызовы. Используйте тот же метод, который вы используете для других вызовов cublas.

  • Вы создаете A в качестве массива векторов. Это не будет работать с cublasSetMatrix. Вместо этого нам нужно создать A как плоский вектор, достаточного размера (N * N), чтобы сохранить всю матрицу.

  • Наконец, cublas ожидает, что матрицы, которые он использует, будут храниться в основном порядке столбца.Если вы передаете массивы C-стиль в строке-мажорный порядке, вы должны использовать транспонирование для этой матрицы в cublasSgemv:

    cublasCheckErrors(cublasSgemv(handle, CUBLAS_OP_T, N, N, &alpha, d_A, N, d_x0, 1, &beta, d_temp, 1)); 
    

Следующий код эти различные проблемы исправлены:

$ cat t235.cu 
#include <cuda.h> 
#include <vector> 
#include <iostream> 
#include <fstream> 
#include <stdio.h> 
#include <stdlib.h> 
#include <cmath> 
#include <cublas_v2.h> 
#include <time.h> 

//#include "timenow.cu" 

// 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) 

// random data filler 
void fillvector(float *data, int N){ 
    for(int i=0; i<N; i++){ 
     data[i] = float(rand() % 10); 
    } 
} 

//printer 
void printer(bool printOut, float *data, int N){ 
    if(printOut == true){ 
    for(int i=0; i<N; i++){ 
     printf("%2.1f ", data[i]); 
    } 
    printf("\n"); 
    } 
} 

///////////////////////////////////////////////////////////////////// 
///////////////////////////////////////////////////////////////////// 

int main(){ 

bool printOut = true; 

int N; 
std::cout << "Enter N: " ; 
std::cin >> N; 

std::vector<float> x0; 
x0.resize(N); 

std::vector<float> p; 
p.resize(N*N); 

// matrix A 
std::vector<float> A(N*N); 
fillvector(A.data(), N*N); 
for (int i=0; i< N; i++){ 
    printer(printOut, &(A[(i*N)]), N); 
    printf("\n");} 
fillvector(x0.data(), N); 
printer(printOut, x0.data(), N); 

printf("\nStarting CUDA computation..."); 
///double startTime = timenow(); 

// device pointers 
float *d_A, *d_x0, *d_temp; 

cudaMalloc((void**)&d_A, N*N*sizeof(float)); 
cudaMalloc((void**)&d_temp, N*sizeof(float)); 
cudaMalloc((void**)&d_x0, N*sizeof(float)); 
cudaCheckErrors("cuda malloc fail"); 

// might need to flatten A... 
cublasCheckErrors(cublasSetVector(N, sizeof(float), &(x0[0]), 1, d_x0, 1)); 
//daMemcpy(d_x0, &x0, N*sizeof(float), cudaMemcpyHostToDevice); 
cublasCheckErrors(cublasSetMatrix(N, N, sizeof(float), &(A[0]), N, d_A, N)); 
//cudaCheckErrors("cuda memcpy of A or x0 fail"); 

float *temp; 
temp = (float *)malloc(N*sizeof(temp)); 

cublasHandle_t handle; 
cublasCheckErrors(cublasCreate(&handle)); 

float alpha = 1.0f; 
float beta = 0.0f; 
cublasCheckErrors(cublasSgemv(handle, CUBLAS_OP_T, N, N, &alpha, d_A, N, d_x0, 1, &beta, d_temp, 1)); 

cublasCheckErrors(cublasGetVector(N, sizeof(float), d_temp, 1, temp, 1)); 
//cudaMemcpy(temp, d_temp, N*sizeof(float), cudaMemcpyDeviceToHost); 
//cudaCheckErrors("returning to host failed"); 

printf("\n"); 
printer(printOut, temp, N); 

/*alpha = -1.0; 
cublasSaxpy(handle, N, &alpha, d_temp, 1, d_v, 1); 
cublasGetVector(N, sizeof(float) * N, d_v, 1, &v, 1); 
printf("\n"); 
for(int i=0; i<N; i++){ 
    printf("%2.1f ",v[i]); 
}*/ 

printf("\nFinished CUDA computations...\n"); 
//double endTime = timenow(); 

//double timeDiff = endTime - startTime; 
//printf("\nRuntime: %2.3f seconds \n", timeDiff); 

cudaFree(d_temp); 
cudaFree(d_A); 
//cudaFree(d_p); 
cudaFree(d_x0); 

return 0; 
} 
$ nvcc -arch=sm_20 -O3 -o t235 t235.cu -lcublas 
$ ./t235 
Enter N: 5 
3.0 6.0 7.0 5.0 3.0 

5.0 6.0 2.0 9.0 1.0 

2.0 7.0 0.0 9.0 3.0 

6.0 0.0 6.0 2.0 6.0 

1.0 8.0 7.0 9.0 2.0 

0.0 2.0 3.0 7.0 5.0 

Starting CUDA computation... 
83.0 86.0 92.0 62.0 110.0 

Finished CUDA computations... 
$ 
+0

Большое спасибо Роберту. Особенно полезны пункты, касающиеся того, как указывать/использовать массивы и их указатели! – rex

+0

Как вы скажете, что это ярмарки против параллельного вычисления процессора, используя, например, 6 ядер? Я спрашиваю, потому что я не совсем уверен, что мой модуль синхронизации работает правильно. Я смотрю N около 1000. – rex

+0

Я также обеспокоен тем, насколько сглаживание матрицы A будет влиять на производительность? есть ли компромисс с точки зрения матричных манипуляций после? – rex