Уменьшите строки матрицы с помощью CUDA
Windows 7, NVidia GeForce 425M.
Я написал простой код CUDA, который вычисляет суммы строк матрицы. Матрица имеет одномерное представление (указатель на число с плавающей точкой).
Серийная версия кода ниже (она имеет 2 петли, как и положено):
void serial_rowSum (float* m, float* output, int nrow, int ncol) {
float sum;
for (int i = 0 ; i < nrow ; i++) {
sum = 0;
for (int j = 0 ; j < ncol ; j++)
sum += m[i*ncol+j];
output[i] = sum;
}
}
Внутри кода CUDA я вызываю функцию ядра, которая подметает матрицу по строкам. Ниже фрагмент вызова ядра:
dim3 threadsPerBlock((unsigned int) nThreadsPerBlock); // has to be multiple of 32
dim3 blocksPerGrid((unsigned int) ceil(nrow/(float) nThreadsPerBlock));
kernel_rowSum<<<blocksPerGrid, threadsPerBlock>>>(d_m, d_output, nrow, ncol);
и функция ядра, которая выполняет параллельную сумму строк (по-прежнему имеет 1 петля):
__global__ void kernel_rowSum(float *m, float *s, int nrow, int ncol) {
int rowIdx = threadIdx.x + blockIdx.x * blockDim.x;
if (rowIdx < nrow) {
float sum=0;
for (int k = 0 ; k < ncol ; k++)
sum+=m[rowIdx*ncol+k];
s[rowIdx] = sum;
}
}
Все идет нормально. Последовательный и параллельный (CUDA) результаты равны.
Все дело в том, что версия CUDA занимает почти вдвое больше времени, чем последовательная, чтобы вычислить, даже если я изменю nThreadsPerBlock параметр: я проверял nThreadsPerBlock от 32 в 1024 (Максимальное количество потоков в блоке разрешено для моей карты).
IMO, размерность матрицы достаточно велика, чтобы оправдать распараллеливание: 90,000 x 1,000,
Ниже я сообщаю, сколько времени прошло для последовательной и параллельной версий, используя разные nThreadsPerBlock, Время сообщается в msec в среднем 100 образцы:
Матрица: nrow = 90000 x ncol = 1000
Серийный номер: среднее время, затраченное на выборку в мсек (100 Образцы): 289.18,
CUDA (32 ThreadsPerBlock): среднее время на выборку в мсек. 100 Образцы): 497.11,
CUDA (1024 ThreadsPerBlock): среднее время на выборку в мсек. 100 Образцы): 699.66,
На всякий случай версия с 32 / 1024nThreadsPerBlock самый быстрый / самый медленный.
Я понимаю, что при копировании с хоста на устройство и наоборот возникают некоторые накладные расходы, но, возможно, медлительность заключается в том, что я не реализую самый быстрый код.
Так как я далеко не эксперт CUDA:
Я кодирую самую быструю версию для этой задачи? Как я могу улучшить свой код? Можно ли избавиться от цикла в функции ядра?
Любые мысли приветствуются.
РЕДАКТИРОВАТЬ 1
Хотя я описываю стандарт rowSum Я заинтересован в AND / OR работа строк, которые имеют (0;1} значения, как rowAND / rowOR, Тем не менее, это не позволяет мне использовать cuBLAS умножить на 1 "s COL Уловка векторного столбца, как предлагают некоторые комментаторы.
РЕДАКТИРОВАТЬ 2
Как подсказывают пользователи других пользователей и здесь одобряются:
Забудьте о попытке написать свои собственные функции, используйте вместо этого библиотеку Thrust, и волшебство придет.
3 ответа
Поскольку вы упомянули, вам нужен общий алгоритм редукции, а не только сумма. Я постараюсь дать 3 подхода здесь. Подход ядра может иметь самую высокую производительность. упорный подход проще всего реализовать. Подход cuBLAS работает только с суммой и имеет хорошие показатели.
Ядро подход
Вот очень хороший документ, рассказывающий, как оптимизировать стандартное параллельное сокращение. Стандартное сокращение можно разделить на 2 этапа.
- Многопоточные блоки каждый уменьшает одну часть данных;
- Один поток блока уменьшается от результата этапа 1 до конечного 1 элемента.
Для вашей задачи мультиредукции (уменьшения рядов мата) достаточно только стадии 1. Идея состоит в том, чтобы уменьшить 1 строку на блок потока. Для дальнейших рассуждений, таких как многострочный блок на одну нить или 1 строку на несколько блоков нити, вы можете обратиться к статье, предоставленной @Novak. Это может улучшить производительность больше, особенно для матриц с плохой формой.
Упорный подход
Общее многократное сокращение может быть сделано thrust::reduction_by_key через несколько минут. Вы можете найти некоторые обсуждения здесь. Определение наименьшего элемента и его положения в каждом столбце матрицы с помощью CUDA Thrust.
тем не мение thrust::reduction_by_key не предполагает, что каждая строка имеет одинаковую длину, поэтому вы получите снижение производительности. Еще один пост Как нормализовать матричные столбцы в CUDA с максимальной производительностью? дает сравнение профилирования между thrust::reduction_by_key и cuBLAS подходят на сумму строк. Это может дать вам общее представление о производительности.
Подход cuBLAS
Сумма строк / столбцов матрицы A может рассматриваться как умножение матрицы на вектор, где все элементы вектора являются единицами. это может быть представлено следующим кодом Matlab.
y = A * ones(size(A,2),1);
где y сумма строк А.
cuBLAS libary предоставляет высокопроизводительную функцию умножения матрицы на вектор cublas<t>gemv() для этой операции.
Результат синхронизации показывает, что эта процедура на 10~50% медленнее, чем простое чтение всех элементов A один раз, что можно рассматривать как теоретический верхний предел производительности для этой операции.
Сокращение строк матрицы может быть решено с помощью CUDA Thrust тремя способами (они могут быть не единственными, но решение этого вопроса выходит за рамки). Как также признается в том же OP, использование CUDA Thrust является предпочтительным для такого рода проблем. Также возможен подход с использованием cuBLAS.
ПОДХОД № 1 -reduce_by_key
Этот подход предлагается на странице примера Thrust. Включает вариант с использованием make_discard_iterator,
ПОДХОД № 2 -transform
Это подход, предложенный Робертом Кровеллой в CUDA Thrust: redu_by_key только для некоторых значений в массиве, основанный на значениях в "ключевом" массиве.
ПОДХОД № 3 -inclusive_scan_by_key
Это подход, предложенный Эриком в Как нормализовать столбцы матрицы в CUDA с максимальной производительностью?,
ПОДХОД № 4 -cublas<t>gemv
Оно использует cuBLASgemv умножить соответствующую матрицу на столбец 1"S.
ПОЛНЫЙ КОД
Вот код, объединяющий два подхода. Utilities.cu а также Utilities.cuh файлы поддерживаются здесь и опускаются здесь. TimingGPU.cu а также TimingGPU.cuh поддерживаются здесь и опущены.
#include <cublas_v2.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/generate.h>
#include <thrust/reduce.h>
#include <thrust/functional.h>
#include <thrust/random.h>
#include <thrust/sequence.h>
#include <stdio.h>
#include <iostream>
#include "Utilities.cuh"
#include "TimingGPU.cuh"
// --- Required for approach #2
__device__ float *vals;
/**************************************************************/
/* CONVERT LINEAR INDEX TO ROW INDEX - NEEDED FOR APPROACH #1 */
/**************************************************************/
template <typename T>
struct linear_index_to_row_index : public thrust::unary_function<T,T> {
T Ncols; // --- Number of columns
__host__ __device__ linear_index_to_row_index(T Ncols) : Ncols(Ncols) {}
__host__ __device__ T operator()(T i) { return i / Ncols; }
};
/******************************************/
/* ROW_REDUCTION - NEEDED FOR APPROACH #2 */
/******************************************/
struct row_reduction {
const int Ncols; // --- Number of columns
row_reduction(int _Ncols) : Ncols(_Ncols) {}
__device__ float operator()(float& x, int& y ) {
float temp = 0.f;
for (int i = 0; i<Ncols; i++)
temp += vals[i + (y*Ncols)];
return temp;
}
};
/**************************/
/* NEEDED FOR APPROACH #3 */
/**************************/
template<typename T>
struct MulC: public thrust::unary_function<T, T>
{
T C;
__host__ __device__ MulC(T c) : C(c) { }
__host__ __device__ T operator()(T x) { return x * C; }
};
/********/
/* MAIN */
/********/
int main()
{
const int Nrows = 5; // --- Number of rows
const int Ncols = 8; // --- Number of columns
// --- Random uniform integer distribution between 10 and 99
thrust::default_random_engine rng;
thrust::uniform_int_distribution<int> dist(10, 99);
// --- Matrix allocation and initialization
thrust::device_vector<float> d_matrix(Nrows * Ncols);
for (size_t i = 0; i < d_matrix.size(); i++) d_matrix[i] = (float)dist(rng);
TimingGPU timerGPU;
/***************/
/* APPROACH #1 */
/***************/
timerGPU.StartCounter();
// --- Allocate space for row sums and indices
thrust::device_vector<float> d_row_sums(Nrows);
thrust::device_vector<int> d_row_indices(Nrows);
// --- Compute row sums by summing values with equal row indices
//thrust::reduce_by_key(thrust::make_transform_iterator(thrust::counting_iterator<int>(0), linear_index_to_row_index<int>(Ncols)),
// thrust::make_transform_iterator(thrust::counting_iterator<int>(0), linear_index_to_row_index<int>(Ncols)) + (Nrows*Ncols),
// d_matrix.begin(),
// d_row_indices.begin(),
// d_row_sums.begin(),
// thrust::equal_to<int>(),
// thrust::plus<float>());
thrust::reduce_by_key(
thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index<int>(Ncols)),
thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index<int>(Ncols)) + (Nrows*Ncols),
d_matrix.begin(),
thrust::make_discard_iterator(),
d_row_sums.begin());
printf("Timing for approach #1 = %f\n", timerGPU.GetCounter());
// --- Print result
for(int i = 0; i < Nrows; i++) {
std::cout << "[ ";
for(int j = 0; j < Ncols; j++)
std::cout << d_matrix[i * Ncols + j] << " ";
std::cout << "] = " << d_row_sums[i] << "\n";
}
/***************/
/* APPROACH #2 */
/***************/
timerGPU.StartCounter();
thrust::device_vector<float> d_row_sums_2(Nrows, 0);
float *s_vals = thrust::raw_pointer_cast(&d_matrix[0]);
gpuErrchk(cudaMemcpyToSymbol(vals, &s_vals, sizeof(float *)));
thrust::transform(d_row_sums_2.begin(), d_row_sums_2.end(), thrust::counting_iterator<int>(0), d_row_sums_2.begin(), row_reduction(Ncols));
printf("Timing for approach #2 = %f\n", timerGPU.GetCounter());
for(int i = 0; i < Nrows; i++) {
std::cout << "[ ";
for(int j = 0; j < Ncols; j++)
std::cout << d_matrix[i * Ncols + j] << " ";
std::cout << "] = " << d_row_sums_2[i] << "\n";
}
/***************/
/* APPROACH #3 */
/***************/
timerGPU.StartCounter();
thrust::device_vector<float> d_row_sums_3(Nrows, 0);
thrust::device_vector<float> d_temp(Nrows * Ncols);
thrust::inclusive_scan_by_key(
thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index<int>(Ncols)),
thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index<int>(Ncols)) + (Nrows*Ncols),
d_matrix.begin(),
d_temp.begin());
thrust::copy(
thrust::make_permutation_iterator(
d_temp.begin() + Ncols - 1,
thrust::make_transform_iterator(thrust::make_counting_iterator(0), MulC<int>(Ncols))),
thrust::make_permutation_iterator(
d_temp.begin() + Ncols - 1,
thrust::make_transform_iterator(thrust::make_counting_iterator(0), MulC<int>(Ncols))) + Nrows,
d_row_sums_3.begin());
printf("Timing for approach #3 = %f\n", timerGPU.GetCounter());
for(int i = 0; i < Nrows; i++) {
std::cout << "[ ";
for(int j = 0; j < Ncols; j++)
std::cout << d_matrix[i * Ncols + j] << " ";
std::cout << "] = " << d_row_sums_3[i] << "\n";
}
/***************/
/* APPROACH #4 */
/***************/
cublasHandle_t handle;
timerGPU.StartCounter();
cublasSafeCall(cublasCreate(&handle));
thrust::device_vector<float> d_row_sums_4(Nrows);
thrust::device_vector<float> d_ones(Ncols, 1.f);
float alpha = 1.f;
float beta = 0.f;
cublasSafeCall(cublasSgemv(handle, CUBLAS_OP_T, Ncols, Nrows, &alpha, thrust::raw_pointer_cast(d_matrix.data()), Ncols,
thrust::raw_pointer_cast(d_ones.data()), 1, &beta, thrust::raw_pointer_cast(d_row_sums_4.data()), 1));
printf("Timing for approach #4 = %f\n", timerGPU.GetCounter());
for(int i = 0; i < Nrows; i++) {
std::cout << "[ ";
for(int j = 0; j < Ncols; j++)
std::cout << d_matrix[i * Ncols + j] << " ";
std::cout << "] = " << d_row_sums_4[i] << "\n";
}
return 0;
}
РЕЗУЛЬТАТЫ ВРЕМЕНИ (проверено на Kepler K20c)
Matrix size #1 #1-v2 #2 #3 #4 #4 (no plan)
100 x 100 0.63 1.00 0.10 0.18 139.4 0.098
1000 x 1000 1.25 1.12 3.25 1.04 101.3 0.12
5000 x 5000 8.38 15.3 16.05 13.8 111.3 1.14
100 x 5000 1.25 1.52 2.92 1.75 101.2 0.40
5000 x 100 1.35 1.99 0.37 1.74 139.2 0.14
Похоже, что подходы № 1 и № 3 превосходят подход № 2, за исключением случаев небольшого количества столбцов. Однако лучшим подходом является подход № 4, который значительно удобнее, чем другие, при условии, что время, необходимое для создания плана, может быть амортизировано во время вычислений.
Если это степень (суммирование строк) операций, которые вам нужно выполнить с этими данными, я бы не ожидал ощутимых преимуществ от графического процессора. У вас есть ровно одна арифметическая операция для каждого элемента данных, и за это вы платите за передачу этого элемента данных в графический процессор. Помимо определенного размера задачи (что бы ни занимало машину), вы не получите дополнительной выгоды от более крупных задач, потому что арифметическая интенсивность равна O (n).
Так что это не особенно захватывающая проблема для GPU.
Но, как показали talonmies, у вас есть проблема слияния в том, как вы ее создали, что еще больше замедлит ход событий. Давайте посмотрим на небольшой пример:
C1 C2 C3 C4
R1 11 12 13 14
R2 21 22 23 24
R3 31 32 33 34
R4 41 42 43 44
Выше приведен простой наглядный пример небольшой части вашей матрицы. Хранение машинных данных таково, что элементы (11), (12), (13) и (14) хранятся в смежных ячейках памяти.
Для объединенного доступа нам нужен такой шаблон доступа, чтобы соседние ячейки памяти запрашивались из одной и той же инструкции, выполняемой через деформацию.
Нам нужно подумать о выполнении вашего кода с точки зрения деформации, то есть 32 потока, выполняющихся в режиме блокировки. Что делает ваш код? Какие элементы он извлекает (запрашивает) на каждом шаге / инструкции? Давайте посмотрим на эту строку кода:
sum+=m[rowIdx*ncol+k];
Смежные потоки в основе имеют смежные (то есть последовательные) значения для rowIdx как вы создали эту переменную. Так когда k = 0, какой элемент данных запрашивается каждым потоком, когда мы пытаемся получить значение m[rowIdx*ncol+k]?
В блоке 0 поток 0 имеет rowIdx из 0. Тема 1 имеет rowIdx 1 и т. д. Таким образом, значения, запрашиваемые каждым потоком в этой инструкции:
Thread: Memory Location: Matrix Element:
0 m[0] (11)
1 m[ncol] (21)
2 m[2*ncol] (31)
3 m[3*ncol] (41)
Но это не коалесцированный доступ! Элементы (11), (21) и т. Д. Не являются смежными в памяти. Для объединенного доступа мы бы хотели, чтобы строка Matrix Element читалась следующим образом:
Thread: Memory Location: Matrix Element:
0 m[?] (11)
1 m[?] (12)
2 m[?] (13)
3 m[?] (14)
Если вы затем работаете в обратном направлении, чтобы определить, что значение ? должно быть, вы придумаете инструкцию примерно так:
sum+=m[k*ncol+rowIdx];
Это даст объединенный доступ, но не даст правильного ответа, потому что мы теперь суммируем столбцы матрицы вместо строк матрицы. Мы можем исправить это, реорганизовав хранилище данных, чтобы они располагались в основном порядке столбцов, а не в порядке строк. (Вы должны быть в состоянии Google, что для идей, верно?) Концептуально, это эквивалентно транспонирования вашей матрицы m, Насколько я понимаю, удобно ли вам это делать или нет - это не входит в сферу ваших вопросов, и на самом деле это не проблема CUDA. Это может быть очень просто для вас, когда вы создаете матрицу на хосте или переносите матрицу с хоста на устройство. Но в целом, я не знаю способа суммировать строки матрицы со 100% объединенным доступом, если матрица хранится в главном порядке строк. (Вы можете прибегнуть к последовательности сокращений строк, но мне это кажется болезненным.)
Когда мы думаем о способах ускорения кода на GPU, нередко стоит подумать о реорганизации нашего хранилища данных для облегчения работы с GPU. Это один из примеров.
И да, то, что я здесь изложил, все еще сохраняет цикл в ядре.
В качестве дополнительного комментария я бы предложил синхронизировать части для копирования данных и части ядра (вычислять) отдельно. По вашему вопросу я не могу сказать, синхронизируете ли вы только ядро или всю (GPU) операцию, включая копии данных. Если вы рассчитываете время копирования данных отдельно, вы можете обнаружить, что только время копирования данных превышает время вашего процессора. Любые усилия по оптимизации вашего кода CUDA не влияют на время копирования данных. Это может быть полезным пунктом данных, прежде чем тратить на это много времени.