4.4 Умножение матриц

Задача. Требуется умножить матрицу на матрицу размерности N * N элементов.

В самой программе необходимо выполнить следующие этапы:

1.
Получить данные для расчетов.
2.
Скопировать эти данные в GPU память.
3.
Произвести вычисление в GPU через функцию ядра.
4.
Скопировать вычисленные данные из GPU памяти в CPU.
5.
Посмотреть результаты.
6.
Высвободить используемые ресурсы.

Переходим непосредственно к написанию кода:

Первым делом напишем функцию ядра, которая и будет осуществлять умножение матриц:

 
// вычисление умножения матрицы на матрицу 
__global__ void mult( int *a, int *b, int *c) { 
    int bx = blockIdx.x; 
    int by = blockIdx.y; 
    int tx = threadIdx.x; 
    int ty = threadIdx.y; 
    float sum = 0.0; 
    int ia=N * 128 * by + N * ty; 
    int ib=128 * bx + tx; 
    for (int k=0; k < N; k++) 
        sum += a[ia + k] * b[ib + k * N]; 
    int ic = N * 128 * by + 128 * bx; 
    c[ic + N * ty + tx] = sum; 
}

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

Чтобы программа представляла рабочий экземплер запишем заголовки и константы.

#include <stdio.h> 
#include <stdlib.h> 
#include <vector> 
#define N (10)

Сразу напишем основные пункты в основном теле программы:

int main( void ) { 
    int *a,*b; 
    int *c; 
    a=(int *)malloc(N*N*sizeof(int)); 
    b=(int *)malloc(N*N*sizeof(int)); 
    c=(int *)malloc(N*N*sizeof(int)); 
    int *dev_a, *dev_b; 
    int *dev_c; 
    cudaMalloc((void**) &dev_a, N*N*sizeof(int)); 
    cudaMalloc((void**) &dev_b, N*N*sizeof(int)); 
    cudaMalloc((void**) &dev_c, N*N*sizeof(int)); 
    for (int i=0; i<N*N; i++) { 
        a[i]=i; 
        b[i]=i*3; 
    } 
    for (int i=0; i<100; i++) c[i]=0; 
    printf("Initialization is finished\n"); 
    cudaMemcpy(dev_a, a,N*N*sizeof(int),cudaMemcpyHostToDevice); 
    cudaMemcpy(dev_b, b,N*N*sizeof(int),cudaMemcpyHostToDevice); 
    cudaMemcpy(dev_c, c,N*N*sizeof(int),cudaMemcpyHostToDevice); 
    mult<<<max(1,N/128),128>>>(dev_a,dev_b,dev_c); 
    printf("Done"); 
    cudaMemcpy(c, dev_c,N*N*sizeof(int),cudaMemcpyDeviceToHost); 
    for (int i=0; i<N*N; i++) { 
        if (i%N==0) printf("\n"); 
        printf("%i ",a[i]); 
    } 
    printf("\n"); 
    for (int i=0; i<N*N; i++) { 
        if (i%N==0) printf("\n"); 
        printf("%i ",b[i]); 
    } 
    printf("\n"); 
    for (int i=0; i<N*N; i++) { 
        if (i%N==0) printf("\n"); 
        printf("%i ",c[i]); 
    } 
    cudaFree( dev_a); 
    cudaFree( dev_b); 
    cudaFree( dev_c); 
    free(a); 
    free(b); 
    free(c); 
    return 0; 
}

В главной части программы мы выполним все основные уже знакомые нам пункты 1-6.