Остаточный расчет с использованием 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 для сокращенного ядра?
Решение
тяга :: 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 а также сопроводительная презентация все еще является хорошим базовым введением для быстрого параллельного сокращения.
Другие решения
Роберт Кровелла уже ответил на этот вопрос и также предложил использовать 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;
}