Остаточный расчет с использованием CUDA

У меня два вектора (oldvector а также newvector). Мне нужно рассчитать значение остатка, который определяется следующим псевдокодом:

residual = 0;
forall i : residual += (oldvector[i] - newvector[i])^2

В настоящее время я рассчитываю это с помощью двух операций CUDA Thrust, которые по существу выполняют:

forall i : oldvector[i] = oldvector[i] - newvector[i]

с последующим thrust::transform_reduce с квадратом в качестве унарного оператора, который делает:

residual = 0;
forall i : residual += oldvector[i]^2;

Проблема с этим, очевидно, является промежуточным хранилищем для глобальной памяти до transform_reduce, Есть ли более эффективный подход к этой проблеме, который бы соединял эти два шага? Помимо написания собственного ядра CUDA, есть ли другой вариант?

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

Если я пишу сокращенное ядро ​​CUDA, были ли сделаны какие-либо улучшения по сравнению с примером CUDA 1.1 для сокращенного ядра?

2

Решение

тяга :: inner_product сделает это за один вызов функции. Ваша оригинальная идея может быть выполнена для работы (сжатие двух векторов и использование thrust::transform_reduce) Этот код демонстрирует оба метода:

#include <iostream>

#include <thrust/tuple.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/transform.h>
#include <thrust/device_vector.h>
#include <thrust/inner_product.h>
#include <thrust/functional.h>

#define N 2

struct zdiffsq{
template <typename Tuple>
__host__ __device__ float operator()(Tuple a)
{
float result = thrust::get<0>(a) - thrust::get<1>(a);
return result*result;
}
};

struct diffsq{
__host__ __device__ float operator()(float a, float b)
{
return (b-a)*(b-a);
}
};

int main(){

thrust::device_vector<float> oldvector(N);
thrust::device_vector<float> newvector(N);
oldvector[0] = 1.0f;  oldvector[1] = 2.0f;
newvector[0] = 2.0f;  newvector[1] = 5.0f;

float result = thrust::inner_product(oldvector.begin(), oldvector.end(), newvector.begin(), 0.0f, thrust::plus<float>(), diffsq());
std::cout << "Result: " << result << std::endl;

float result2 = thrust::transform_reduce(thrust::make_zip_iterator(thrust::make_tuple(oldvector.begin(), newvector.begin())), thrust::make_zip_iterator(thrust::make_tuple(oldvector.end(), newvector.end())), zdiffsq(), 0.0f, thrust::plus<float>());
std::cout << "Result2: " << result2 << std::endl;
}

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

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

И да, Параллельная редукция CUDA а также сопроводительная презентация все еще является хорошим базовым введением для быстрого параллельного сокращения.

3

Другие решения

Роберт Кровелла уже ответил на этот вопрос и также предложил использовать CUB.

В отличие от Thrust, CUB оставляет критичные для производительности параметры, такие как выбор конкретного алгоритма сокращения, который будет использоваться, и степень параллелизма, не зависящая от пользователя, выбираемая пользователем. Эти параметры можно настроить, чтобы максимизировать производительность для конкретной архитектуры и приложения. Параметры могут быть указаны во время компиляции, что позволяет избежать снижения производительности во время выполнения.

Ниже приведен полностью проработанный пример того, как использовать CUB для расчета остатков.

#include <cub/cub.cuh>
#include <cuda.h>

#include "Utilities.cuh"
#include <iostream>

#define BLOCKSIZE           256
#define ITEMS_PER_THREAD    8

const int N = 4096;

/******************************/
/* TRANSFORM REDUCTION KERNEL */
/******************************/
__global__ void TransformSumKernel(const float * __restrict__ indata1, const float * __restrict__ indata2, float * __restrict__ outdata) {

unsigned int tid = threadIdx.x + blockIdx.x * gridDim.x;

// --- Specialize BlockReduce for type float.
typedef cub::BlockReduce<float, BLOCKSIZE * ITEMS_PER_THREAD> BlockReduceT;

__shared__ typename BlockReduceT::TempStorage  temp_storage;

float result;
if(tid < N) result = BlockReduceT(temp_storage).Sum((indata1[tid] - indata2[tid]) * (indata1[tid] - indata2[tid]));

if(threadIdx.x == 0) atomicAdd(outdata, result);

return;
}

/********/
/* MAIN */
/********/
int main() {

// --- Allocate host side space for
float *h_data1      = (float *)malloc(N * sizeof(float));
float *h_data2      = (float *)malloc(N * sizeof(float));
float *h_result     = (float *)malloc(sizeof(float));

float *d_data1;     gpuErrchk(cudaMalloc(&d_data1, N * sizeof(float)));
float *d_data2;     gpuErrchk(cudaMalloc(&d_data2, N * sizeof(float)));
float *d_result;    gpuErrchk(cudaMalloc(&d_result, sizeof(float)));

for (int i = 0; i < N; i++) {
h_data1[i] = 1.f;
h_data2[i] = 3.f;
}

gpuErrchk(cudaMemcpy(d_data1, h_data1, N * sizeof(float), cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_data2, h_data2, N * sizeof(float), cudaMemcpyHostToDevice));

TransformSumKernel<<<iDivUp(N, BLOCKSIZE), BLOCKSIZE>>>(d_data1, d_data2, d_result);
gpuErrchk(cudaPeekAtLastError());
gpuErrchk(cudaDeviceSynchronize());

gpuErrchk(cudaMemcpy(h_result, d_result, sizeof(float), cudaMemcpyDeviceToHost));

std::cout << "output: ";
std::cout << h_result[0];
std::cout << std::endl;

gpuErrchk(cudaFree(d_data1));
gpuErrchk(cudaFree(d_data2));
gpuErrchk(cudaFree(d_result));

return 0;
}
1