4.3 Скалярное произведение векторов

Задача. Требуется вычислить скалярное произведение векторов размерностью N элементов.

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

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

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

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

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

// процедура вычисления скалярного умножения на видеокарте, возвращает массив по количеству блоков с суммой внутри блока 
__global__ void dot( float *a, float *b, float *c ) { 
    __shared__ float cache[threadsPerBlock]; 
    int tid = threadIdx.x + blockIdx.x * blockDim.x; 
    int cacheIndex = threadIdx.x; 
    float temp = 0; 
    while (tid < N) { 
        temp += a[tid] * b[tid]; 
        tid += blockDim.x * gridDim.x; 
    } 
    // сохраняем в кеш 
    cache[cacheIndex] = temp; 
    // ждем выполнения вычислений на всех нитях данного блока 
    __syncthreads(); 
    // редукция /2 
    int i = blockDim.x / 2; 
    while (i != 0) { 
        if (cacheIndex < i) cache[cacheIndex] += cache[cacheIndex + i]; 
        __syncthreads(); 
        i /= 2; 
    } 
    if (cacheIndex == 0) c[blockIdx.x] = cache[0]; 
}

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

#include <stdio.h> 
#define imin(a,b) (a<b?a:b) 
 
const int N = 60 * 1024; 
const int threadsPerBlock = 256; 
const int blocksPerGrid = imin( 32, ( (N+threadsPerBlock - 1) / threadsPerBlock) );

Выделим отдельную функцию на хосте, в которой будут инициализироваться данные, выделяться память на видеокарте, копироваться данные в данные на видеокарте, и собственно запускаться функция ядра:

 
// функция суммирования по блокам на центральном процессоре 
 
float vectorbyvector(float *dev_a, float *dev_b, float *dev_c) { 
   float *partial_c; 
   float *dev_a, *dev_b; 
   partial_c = (float *)malloc(blocksPerGrid*sizeof(float)); 
   dot<<<blocksPerGrid,threadsPerBlock>>>(dev_a,dev_b,dev_c); 
   cudaMemcpy( partial_c, dev_c, blocksPerGrid*sizeof(float), cudaMemcpyDeviceToHost); 
   float dotres = 0; 
   for (int i=0; i<blocksPerGrid; i++) dotres+=partial_c[i]; 
   delete [] partial_c; 
   return dotres; 
}

Данная функция выполняет пункты 2, 3, 4 и 6. В начале функции происходит вычисление количества блоков, который зависит от мощности вектора и количества элементов в блоке. Посредством cudaMalloc выделяем память в видеокарте. Переносим данные в память видеокарты с помошью cudaMemcpy. Запускаем вычисления на видеокарте с помошью вызова функции ядра summOfVector. Заметим, что следует сделать синхронизацию нитей. После вычислений на видеокарте, остаеться просуммировать элементы куда более меньшего массива на хосте. В конце функции происходит освобождение использованной памяти. И в качестве результата данная функция возвращает скаляр, который является суммой всех элементов вектора.

int main( void ) { 
    printf("Initialization... "); 
    int *inputdata,outputdata; 
    inputdata=(int *)malloc(N*sizeof(int)); 
    for (int i=0; i<N; i++) { 
        inputdata[i]=i; 
    } 
    printf("Done!\n"); 
    outputdata=summOfVector(inputdata); 
    printf("Result: %i",outputdata); 
    free(inputdata); 
    return 0; 
}

В главной части программы остаеться лишь инициализировать данные, запустить предыдущую функцию и вывести результаты.