2014-09-28 2 views
0

У меня есть 200 матриц A [i] (размерность 4096 * 48) и 48 векторов v [j] (размерность 48 * 1). Я хочу рассчитать A [i] * v [j], (i = 0: 199, j = 1: 47).Размер сетки и размер блока

Я думаю о том, как упорядочить размер сетки и размер блока со вчерашнего дня. Но сейчас я не понимаю ответа. Может ли кто-нибудь дать мне совет?

Максимальное количество блоков - 512. Это моя рабочая среда. enter image description here

Следующий мой код. Он работает правильно. Я проверил. Но это медленнее, чем Matlab :(

#include<iostream> 
#include <mat.h> 
#include <time.h> 
#include <cuda_runtime.h> 
#include "cuda.h" 

using std::cout; 
using std::endl; 
using namespace cv; 
using namespace std; 

#include <limits> 
#include <iostream> 
#include <cstdlib> 
using namespace std; 

#define kernel_size 48 

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

typedef struct { 
    int width; 
    int height; 
    int stride; 
    float* elements; 
} Matrix; 



// Forward declaration of the matrix multiplication kernel 
__global__ void MatMulKernel(const Matrix, const Matrix, Matrix); 
// Matrix multiplication - Host code 
// Matrix dimensions are assumed to be multiples of BLOCK_SIZE 
void MatMul(const Matrix A, const Matrix B, Matrix C) 
{ 
    // Load A and B to device memory 
    Matrix d_A; 
    d_A.width = d_A.stride = A.width; d_A.height = A.height; 
    size_t size = A.width * A.height * sizeof(float); 
    cudaMalloc(&d_A.elements, size); 
    cudaMemcpy(d_A.elements, A.elements, size, 
     cudaMemcpyHostToDevice); 
    Matrix d_B; 
    d_B.width = d_B.stride = B.width; d_B.height = B.height; 
    size = B.width * B.height * sizeof(float); 
    cudaMalloc(&d_B.elements, size); 
    cudaMemcpy(d_B.elements, B.elements, size, 
     cudaMemcpyHostToDevice); 
    // Allocate C in device memory 
    Matrix d_C; 
    d_C.width = d_C.stride = C.width; d_C.height = C.height; 
    size = C.width * C.height * sizeof(float); 
    cudaMalloc(&d_C.elements, size); 
    // Invoke kernel 
    dim3 dimBlock(1,B.height); 
    dim3 dimGrid(A.height, C.width); 
    MatMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C); 
    // Read C from device memory 
    cudaMemcpy(C.elements, d_C.elements, size, 
     cudaMemcpyDeviceToHost); 
    // Free device memory 
    cudaFree(d_A.elements); 
    cudaFree(d_B.elements); 
    cudaFree(d_C.elements); 
} 
// Matrix multiplication kernel called by MatMul() 
__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C) 
{ 
    // Block row and column 
    int blockCol = blockIdx.y; 
    int blockRow = blockIdx.x; 

    float Cvalue = 0; 
    // Thread row and column within Csub 
    int row = threadIdx.y; 
    int col = threadIdx.x; 
    // Loop over all the sub-matrices of A and B that are 
    // required to compute Csub 
    // Multiply each pair of sub-matrices together 
    // and accumulate the results 


    // Shared memory used to store Asub and Bsub respectively 
    __shared__ float As[1][kernel_size]; 
    __shared__ float Bs[kernel_size][1]; 
    // Load Asub and Bsub from device memory to shared memory 
    // Each thread loads one element of each sub-matrix 


    As[0][row] = A.elements[blockRow * A.stride + row+B.height*blockCol]; 
    Bs[row][0] = B.elements[row]; 
    // Synchronize to make sure the sub-matrices are loaded 
    // before starting the computation 
    __syncthreads(); 
    // Multiply Asub and Bsub together 
    for (int e = 0; e < B.height; ++e) 
    { 
     Cvalue += As[0][e] * Bs[e][0]; 

    } 
    // Synchronize to make sure that the preceding 
    // computation is done before loading two new 
    // sub-matrices of A and B in the next iteration 
    __syncthreads(); 

    // Write Csub to device memory 
    // Each thread writes one element 
    C.elements[blockRow * C.stride +blockCol]= Cvalue; 
} 

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





float * gen_matrix(int n /*row*/, int m /*col*/){ 

    float *A; 
    //srand(1023); 
    A = (float *) malloc(n*m*sizeof(float)); 

    for(int row = 0;row < n;row++) 
     for(int col = 0;col < m;col++) { 
      A[row*m+col] = rand()%10; 
     } 

     /* 
     // print matrix elements. 
     for (int i = 0; i < n; ++i) { 
     for (int j = 0; j < m; ++j) 
     cout << " [" << i << "," << j << "] " << A[i*m+j] ; 
     cout << endl; 
     } 
*/ 
     return A; 
} 



int main() 
{ 
    int k=kernel_size; 
    int s=2000; 
    int m =4096; 
    //int m=2; 
    //int s=1; 
    int n = k*s; 
    float *Ae = gen_matrix(m,n); 
    float *Be= gen_matrix(k,1);00 
    float *Ce=(float *) malloc(m*s*sizeof(float)); 

    Matrix A ={n,m,n,Ae}; 
    Matrix B ={1,k,1,Be}; 
    Matrix C ={s,m,s,Ce}; 

    const clock_t begin_time = clock(); 
    MatMul(A, B, C); 
    std::cout << float(clock() - begin_time)/CLOCKS_PER_SEC; 

    for (int i = 0; i < 3; ++i) { 
     for (int j = 0; j <7; ++j) 
      cout << " [" << i << "," << j << "] " << Ce[i*m+j] ; 
     cout << endl; 
    } 


    //check 
    float *Ce2=(float *) malloc(s*m*sizeof(float)); 
    for (int i = 0; i < m; i++) 
    { 
     for (int j = 0; j < s; j++) 
     { 
      Ce2[i*s+j]=0; 
     } 
    } 
    for (int i = 0; i < m; i++) 
    { 
     for (int j = 0; j < s; j++) 
     { 
      for (int ind = 0; ind < k; ind++) 
      { 
       Ce2[i*s+j]=Ce2[i*s+j]+Ae[j*k+ind+i*k*s]*Be[ind]; 
      // printf("%f---****%f\n",Ae[j*k+ind+i*k*s],Be[ind]); 
      } 
      if (Ce2[i*s+j]!= Ce[i*s+j]) 
      { 
       printf("%f----%f\n",Ce2[i*s+j],Ce[i*s+j]); 
      } 

     } 

    } 





    free(Ae); 
    free(Be); 
    free(Ce); 
} 
+1

Забудьте о своей структуре Matrix и подумайте, как вы можете упорядочить свои данные в одномерных массивах. Как только вы это сделаете, размеры сетки и блока в этом случае кажутся почти произвольными. –

ответ

1

Это просто проблема умножения матрицы на матрицу. Если вы хотите вещи, чтобы быстро бегать, вы не должны писать свой собственный код матрицы-матрицы умножения. Использование CUBLAS Sgemm.

Концептуально, если вы устраиваете свои A матрицы, как это:

[A0] 
[A1] 
[A2] 
... 
[A199] 

тогда вы будете иметь новую матрицу AA, которая (4096 * 200) строк х 48 столбцов

.

Упорядочивание 48 V векторов (48x1) в матрице 48х48 (VV):

[V0][V1][V2]...[V47] 

(каждый V вектор представляет собой столбец новой матрицы VV)

теперь у вас есть один матрица проблема умножения (AA * VV), то есть (4096 * 200) x48, умноженная на 48x48, что дает результат (4096 * 200) x 48. Этот результат имеет один вектор-столбец длиной 4096 * 200, который содержит 200 результатов отдельных умножений матричных векторов, которые вы пытались сделать. 200 результатов на столбец * 48 столбцов объединяются, чтобы дать вам все результаты, которые создала ваша оригинальная проблема. Первый столбец будет содержать результаты [V0] умноженных на каждом из 200 A матриц, второй столбец будет содержать результаты [V1] умноженных на каждом из 200 A матриц и т.д.

После того, как вы расположили свои данные, как это , использование CUBLAS Sgemm должно быть самым быстрым подходом на GPU. Обратите внимание, что CUBLAS ожидает, что базовое хранилище будет иметь значение столбца, поэтому, если вы переставляете свои данные, вы, вероятно, захотите это учитывать. Существует CUDA sample code for CUBLAS matrix multiplication.

В вашем коде у вас на самом деле есть 2000 A матриц, но ваш вопрос относится к 200. Я использовал 200, например, в своем ответе, но концепция будет такой же, как и с матрицами 2000 A.

Смежные вопросы