4.2 Сумма элементов вектора

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

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

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

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

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

// Функция сложения элементов вектора 
 
__global__ void summOfVector(int *a, int *b) { 
    __shared__ int data[BLOCK_SIZE]; // создание массива в Shared Memory девайса 
    int tid=threadIdx.x; 
    int idx=blockIdx.x*blockDim.x+threadIdx.x; 
    data[tid]=(idx<N) ? a[idx] : 0; // копируем с глобальной памяти 
    //данные в Shared Memory, вызов которой намного быстрее. 
    __syncthreads();  // ждем пока все нити(потоки) скопируют данные. 
 
    for(int s = blockDim.x/2; s>0; s=s/2) { 
        if (tid<s) data[tid]+=data[tid+s]; 
        __syncthreads(); 
    } //суммируем так, чтобы нити(потоки) не работали с одной и той же памятью. 
    if (tid==0) 
        b[blockIdx.x]=data[0]; 
    // после прохода цикла, сумма всех элементов вектора a, которое вычисляется 
    // в данном блоке будет находится в data[0] 
    // по этому сохраняем его в глобальной памяти. 
}

Таким образом, мы будем использовать более быструю разделяемую память.Как нам известно разделяемая память является памятью общего доступа для нитей внутри блока. В этой процедуре мы применяем распространенный прием редукцию.Редукция – это процедура, в результате которой из входного массива получается массив меньшей размерности. Идея в том, что каждая нить складывает два значения находящиеся в массиве data[], и помещает результат снова в data[]. Поскольку при этом два элемента объединяются в один, то по завершении данного шага количество элементов уменьшится вдвое. На следующем шаге та же операция применяется к оставшейся половине.

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

#include <stdio.h> 
#include <stdlib.h> 
#define N (10000) 
#define BLOCK_SIZE 512

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

// инициализация векторов 
void initvector() { 
    vecTemp = (float*)malloc(N*sizeof(float)); 
    vecX = (float*)malloc(N * sizeof(float)); 
    vecP = (float*)malloc(N * sizeof(float)); 
    vecR = (float*)malloc(N * sizeof(float)); 
    vecB = (float*)malloc(N * sizeof(float)); 
 
    for (int i = 0; i < N; i++) { 
        vecX[i] = vecP[i] = vecR[i] = vecTemp[i] = 0.0; 
        vecB[i] = 1.0 / N; 
        vecP[i] = 1.0; 
    } 
}

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

 
float floatvectorbyvector(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); 
 
floatdotres = 0; 
for (int i=0; i<blocksPerGrid; i++) dotres += partial_c[i]; 
 
delete [] partial_c; 
 
return dotres; 
}

Данная функция выполняет роль интегратора. И в качестве результата данная функция возвращает скаляр, который является суммой всех элементов выходного вектора.

int main( void ) { 
    initvector(); 
 
    float *dev_a, *dev_b, *dev_c; 
 
    cudaMalloc( (void**)&dev_a, N*sizeof(float) ); 
    cudaMalloc( (void**)&dev_b, N*sizeof(float) ); 
    cudaMalloc( (void**)&dev_c, blocksPerGrid*sizeof(float)); 
 
    cudaMemcpy( dev_a, vecB, N * sizeof(float), cudaMemcpyHostToDevice); 
    cudaMemcpy( dev_b, vecP, N * sizeof(float), cudaMemcpyHostToDevice); 
 
    alpha = vectorbyvector(dev_a, dev_b, dev_c); 
 
    cudaFree( dev_a); 
    cudaFree( dev_b); 
    cudaFree( dev_c); 
 
    printf("Результат скалярного умножения %f \n", alpha); 
}

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