2013-08-15 2 views
0

То, что я пытаюсь сделать, это Multiply Matrix A & Матрица B, а затем из матрицы продукта. Я получаю индекс максимального значения за столбец. Но, к сожалению, только первые 128 * 128 значений матричного умножения являются правильными, а другие - просто мусором. Я не совсем понимаю, как это работает. Я прошу вас любезно провести меня с этим.Матрица Умножение, дающее неправильный результат

#include<stdio.h> 
#include "cuda.h" 
#include<stdlib.h> 

#define blockD 32 
const int wA = 128; 
const int hA = 4096;  
const int wB = 4096; 
const int hB = wA; 

main(void){ 

    void MatrixMultiplication(float *, float *, float *, float *); 

    int size_A = wA * hA * sizeof(float); 
    int size_B = wB * hB * sizeof(float); 
    int size_C = wB * hA * sizeof(float); 
    int size_max = 2 * wB * sizeof(float); 
    float *M, *N, *P, *C; 

    // allocate memory on the CPU 
    M = (float*)malloc(size_A); 
    N = (float*)malloc(size_B); 
    P = (float*)malloc(size_max); 
    C = (float*)malloc(size_C); 

    // initialize the matrices 
    for (int y=0; y < hA; y++) { 
     for (int x=0; x < wA; x++){ 
      M[y*wA + x] = 32; //x + y*wA; 
     } 
    } 

    for (int y=0; y<hB; y++) { 
     for (int x=0; x<wB; x++){ 
      N[y*wB + x] = 21; //x + y*wB; 
     } 
    } 


    MatrixMultiplication(M, N, P, C); 

    //Write 
    FILE *f1; 
    int i,j; 
    f1 = fopen("C.txt","w"); 
    for(i = hA - 2 ; i < hA; i ++){ 
    for(j = 0; j < wB; j++){ 
     fprintf(f1,"%d\t",int(C[i*wB + j])); 
    } 
    fprintf(f1,"\n"); 
    } 
    fclose(f1); 

    // free the memory allocated on the CPU 
    free(M); 
    free(N); 
    free(P); 
    free(C); 
    cudaDeviceReset(); 
    return 0; 
} 


__device__ void MaxFunction(float* Pd, float* max) 
{ 
int x = (threadIdx.x + blockIdx.x * blockDim.x); 
int y = (threadIdx.y + blockIdx.y * blockDim.y); 

int k = 0; 

int temp = 0; int temp_idx = 0; 
for (k = 0; k < wB; ++k) { 
      if(Pd[x*wB + k] > temp){ 
       temp = Pd[x*wB + k]; 
       temp_idx = x*wB + k; 
      } 
    } 
    max[y*2 + 0] = temp; 
    max[y*2 + 1] = temp_idx; 
} 


__global__ void MatrixMulKernel(float* Md, float* Nd, float* Pd, float* max) 
{ 
    // declare cache in the shared memory 
    __shared__ float Mds[blockD][blockD]; 
    __shared__ float Nds[blockD][blockD]; 

    float Pvalue = 0; 
    // Loop over the Md and Nd block dimension required to compute the Pd element 
    for (int m = (wA * blockD * blockIdx.y), n = (blockD * blockIdx.x); 
          m < ((wA * blockD * blockIdx.y)+wA-1); 
             m += blockD, n += (blockD*hB)){ 

    // collaboratively loading of Md and Nd blocks into shared memory  
    Mds[threadIdx.y][threadIdx.x] = Md[m + wA * threadIdx.y + threadIdx.x]; 
    Nds[threadIdx.y][threadIdx.x] = Nd[n + wA * threadIdx.y + threadIdx.x]; 
    __syncthreads(); 

    // keep track of the running sum  
    for (int k = 0; k < blockD; k++) 
     Pvalue += Mds[threadIdx.y][k] * Nds[k][threadIdx.x]; 
    __syncthreads(); 
    } 

    // write back to the global memory 
    int p = hB * blockD * blockIdx.y + blockD * blockIdx.x; 
    Pd[p + hB * threadIdx.y + threadIdx.x] = Pvalue; 
    __syncthreads(); 

    MaxFunction(Pd, max); 

} 

void MatrixMultiplication(float *M, float *N, float *P, float *C) { 

    int size_A = wA * hA * sizeof(float); 
    int size_B = wB * hB * sizeof(float); 
    int size_C = wB * hA * sizeof(float); 
    int size_max = 2 * wB * sizeof(float); 
    float *Md, *Nd, *Pd, *max; 

    // allocate memory on the GPU 
    cudaMalloc((void**)&Md, size_A); 
    cudaMalloc((void**)&Nd, size_B); 
    cudaMalloc((void**)&Pd, size_C); 
    cudaMalloc((void**)&max, size_max); 

    // transfer M and N to device memory 
    cudaMemcpy(Md, M, size_A, cudaMemcpyHostToDevice); 
    cudaMemcpy(Nd, N, size_B, cudaMemcpyHostToDevice); 

    // kernel invocation code 
    dim3 dimBlock(blockD, blockD); 
    dim3 dimGrid(wA/blockD, hB/blockD); 

    //Execute Kernel 
    MatrixMulKernel<<<dimGrid, dimBlock>>>(Md, Nd, Pd, max); 

    // transfer P from device  
    cudaMemcpy(P, max, size_max, cudaMemcpyDeviceToHost); 
    cudaMemcpy(C, Pd, size_C, cudaMemcpyDeviceToHost); 

    // free the memory allocated on the GPU 
    cudaFree(Md); 
    cudaFree(Nd); 
    cudaFree(Pd); 
    cudaFree(max); 
} 
+0

Это точно такой же код и тот же вопрос, что и в предыдущем вопросе. Пожалуйста, не повторяйте тот же вопрос еще раз. – talonmies

+0

Я согласен с тем же кодом. Но я не могу найти ответ. – Krsrb1

+0

Это не повод для публикации дублирующего вопроса. Ключом к получению помощи является редактирование существующего вопроса, чтобы было легче ответить. Прямо сейчас у вашего кода есть две отдельные проблемы - матричное умножение и сокращение. Выберите проблему. Улучшите код - например, я не вижу проверки ошибок API CUDA. Вы даже уверены, что код действительно работает до завершения? Используйте предоставленные инструменты - отладчик, cuda-memcheck. Улучшите вопрос с тем, что вы найдете - [SO] не является бесплатной службой отладки, где мы делаем вашу работу за вас. Помогите нам помочь вам ... – talonmies

ответ

1

В вашем коде у вас, кажется, есть несколько проблем. Одна из проблем заключается в том, вместо этого:

dim3 dimGrid(wA/blockD, hB/blockD); 

Вы должны иметь это:

dim3 dimGrid(wB/blockD, hA/blockD); 

В конечном итоге вы должны один поток в сетке для каждой точки выхода. Ваша формулировка давала вам сетку из 4 блоков на 4 блока, тогда как вам нужна сетка из 128 блоков на 128 блоков.

Другая проблема, которую я нашел с кодом в этих строках в ядре:

int p = hB * blockD * blockIdx.y + blockD * blockIdx.x; 
Pd[p + hB * threadIdx.y + threadIdx.x] = Pvalue; 

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

Pd[(threadIdx.x + (blockIdx.x * blockDim.x)) + ((threadIdx.y + (blockIdx.y * blockDim.y))*(gridDim.x*blockDim.x))] = Pvalue; 

Когда я сделал две вышеупомянутые изменения в свой код, я получил то, что я считаю, это правильные результаты во всем массиве. И потребовалось около 32 секунд на моей машине, чтобы запустить его. (Обратите внимание, что я не пробовал исправить ваш исходный код максимального поиска - см. Ниже, чтобы получить лучший подход.)

Основываясь на вашем предыдущем вопросе, вы, похоже, беспокоились о скорости. Если вы хотите быстро умножить матрицу, вы должны использовать cublas. Следующий код показывает, как использовать cublas для умножения двух обычных C-образных матриц (они не обязательно должны быть квадратными). Я также включил ядро ​​поиска столбца-max, которое будет быстрым, когда количество столбцов велико (скажем, более 500 или около того. В вашем примере есть 4096 столбцов). Для небольшого числа столбцов могут быть более быстрые способы выполнения этой функции, но небольшое количество столбцов также указывает на то, что общий размер проблемы может быть небольшим, и поэтому скорость (этого фрагмента кода) на самом деле не будет проблемой.

Вот код:

#include <stdio.h> 
#include <cublas_v2.h> 
#define VERBOSE 1 
#define nTPB 64 
#define ROW_A 4 
#define COL_A 4 
#define ROW_B COL_A 
#define COL_B 4 
#define ROW_C ROW_A 
#define COL_C COL_B 
#define SIZ_A (ROW_A*COL_A) 
#define SIZ_B (ROW_B*COL_B) 
#define SIZ_C (ROW_C*COL_C) 



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

__global__ void col_max(float *mat, float *max, unsigned int *midx, unsigned int rows, unsigned int cols){ 
    int idx = threadIdx.x + blockDim.x*blockIdx.x; 
    if (idx < cols){ 
    float tempmax = mat[idx]; 
    unsigned int tempmidx = 0; 
    for (int i = 1; i< rows; i++) 
     if (mat[idx + (i*cols)] > tempmax){ 
     tempmax = mat[idx + (i*cols)]; 
     tempmidx = i;} 
    max[idx] = tempmax; 
    midx[idx] = tempmidx; 
    } 
} 

int main(){ 

    float *h_A, *h_B, *h_C, *d_A, *d_B, *d_C, *h_max, *d_max; 
    unsigned int *h_idx, *d_idx; 

    h_A = (float *)malloc(SIZ_A*sizeof(float)); 
    if (h_A==0) {printf("malloc fail\n"); return -1;} 
    h_B = (float *)malloc(SIZ_B*sizeof(float)); 
    if (h_B==0) {printf("malloc fail\n"); return -1;} 
    h_C = (float *)malloc(SIZ_C*sizeof(float)); 
    if (h_C==0) {printf("malloc fail\n"); return -1;} 
    h_max = (float *)malloc(COL_C*sizeof(float)); 
    if (h_max==0) {printf("malloc fail\n"); return -1;} 
    h_idx = (unsigned int*)malloc(COL_C*sizeof(unsigned int)); 

    if (h_idx==0) {printf("malloc fail\n"); return -1;} 

    cudaMalloc((void **)&d_A, SIZ_A*sizeof(float)); 
    cudaMalloc((void **)&d_B, SIZ_B*sizeof(float)); 
    cudaMalloc((void **)&d_C, SIZ_C*sizeof(float)); 
    cudaMalloc((void **)&d_max, COL_C*sizeof(float)); 
    cudaMalloc((void **)&d_idx, COL_C*sizeof(unsigned int)); 
    cudaCheckErrors("cuda malloc fail"); 

    // initialize data 
    for (int i=0; i< SIZ_A; i++) h_A[i] = (float)(i+1); 
    for (int i=0; i< SIZ_B; i++) h_B[i] = (float)(i+2); 

    cudaMemcpy(d_A, h_A, SIZ_A*sizeof(float), cudaMemcpyHostToDevice); 
    cudaMemcpy(d_B, h_B, SIZ_B*sizeof(float), cudaMemcpyHostToDevice); 
    cudaCheckErrors("cuda memcpy 1 fail"); 
    const float alpha = 1.0f; 
    const float beta = 0.0f; 
    cublasHandle_t handle; 
    cublasCheckErrors(cublasCreate(&handle)); 
    // C = A*B 
    // due to cublas expecting column-major storage, parameters 
    // are scrambled 
    cublasCheckErrors(cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, COL_B, ROW_A, COL_A, &alpha, d_B, COL_B, d_A, COL_A, &beta, d_C, COL_C)); 
    cudaMemcpy(h_C, d_C, SIZ_C*sizeof(float), cudaMemcpyDeviceToHost); 
    cudaCheckErrors("cuda memcpy 2 fail"); 
    col_max<<<(COL_C + nTPB - 1)/nTPB, nTPB>>>(d_C, d_max, d_idx, ROW_C, COL_C); 
    cudaCheckErrors("kernel launch fail"); 
    cudaMemcpy(h_max, d_max, COL_C*sizeof(float), cudaMemcpyDeviceToHost); 
    cudaMemcpy(h_idx, d_idx, COL_C*sizeof(unsigned int), cudaMemcpyDeviceToHost); 
    cudaCheckErrors("cuda memcpy 3 fail/kernel fail"); 

    if (VERBOSE){ 
    printf("A: \n"); 
    for (int i=0; i< ROW_A; i++){ 
     for (int j=0; j< COL_A; j++) 
     printf("%7.5G", h_A[j+(i*COL_A)]); 
     printf("\n");} 
    printf("B: \n"); 
    for (int i=0; i< ROW_B; i++){ 
     for (int j=0; j< COL_B; j++) 
     printf("%7.5G", h_B[j+(i*COL_B)]); 
     printf("\n");} 
    printf("C = A*B: \n"); 
    for (int i=0; i< ROW_C; i++){ 
     for (int j=0; j< COL_C; j++) 
     printf("%7.5G", h_C[j+(i*COL_C)]); 
     printf("\n");} 
    printf("COLUMN MAX:\n"); 
    for (int i=0; i< COL_C; i++) 
     printf("%7.5G", h_max[i]); 
    printf("\nCOLUMN MAX IDX:\n"); 
    for (int i=0; i< COL_C; i++) 
     printf("%7d", h_idx[i]); 
    } 
    printf("\n finished!\n"); 
    return 0; 
} 

Вот что я использовал для компиляции:

$ nvcc -arch=sm_20 -O3 -o t221 t221.cu -lcublas 

А вот пример вывода:

$ cuda-memcheck ./t221 
========= CUDA-MEMCHECK 
A: 
     1  2  3  4 
     5  6  7  8 
     9  10  11  12 
    13  14  15  16 
B: 
     2  3  4  5 
     6  7  8  9 
    10  11  12  13 
    14  15  16  17 
C = A*B: 
    100 110 120 130 
    228 254 280 306 
    356 398 440 482 
    484 542 600 658 
COLUMN MAX: 
    484 542 600 658 
COLUMN MAX IDX: 
     3  3  3  3 
finished! 
========= ERROR SUMMARY: 0 errors 
$ 

Когда я протянул код для обработки те же размеры, которые вы указали, (A = 4096x128, B = 128x4096), мне потребовалось около 1 секунды на моей машине. Так что это намного быстрее, чем ваш код. Однако, когда я беру ваш код и прокомментирую ваш вызов MaxFunction в ядре, для вычисления результата умножения матрицы потребуется всего 1 секунду. Поэтому, если вы хотите сохранить код умножения на матрицу (т. Е. Не использовать cublas), вы можете разбить код на 2 ядра и использовать свою процедуру умножения в первом ядре с моей программой максимального поиска (col_max) во втором ядре, а также вероятно, получить довольно быстрый результат.

Как указано в @talonmies, если вы работаете на машине с Windows, убедитесь, что знаете о ветвлениях TDR Windows.(поиск в окне поиска в правом верхнем углу, если необходимо)

+0

Стоит отметить, что код умножения матрицы, размещенный в исходном вопросе, действительно работает нормально. Я подозреваю, что он запускается на медленном устройстве и нажимает сторожевой таймер дисплея. На самом деле здесь нет никаких вопросов, но спасибо за сообщение разумного ответа в любом случае ... – talonmies

+0

Теперь я отредактировал свой ответ с моими исправлениями, чтобы получить код OP, размещенный в этом вопросе, чтобы генерировать (я думаю) правильную матрицу, умножать результаты. Я довольно убежден, что код OP в этом вопросе не дает правильных результатов умножения матрицы. –

+0

Пока матрицы являются квадратными (wA = wB = hB) и круглыми кратными размеру плитки (так 32), код умножения матрицы работал при каждом размере, в котором я пробовал его от 128 до 4096. Легко проверить, каждая запись должна быть wA * 32 * 21. Это повторяется снова и снова, код умножения матрицы SDK получает неправильное использование, а затем возникают вопросы/жалобы о том, почему он не работает ... – talonmies

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