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(32ThreadsPerBlock):毫秒(每个样品经过的平均时间100样品)497.11.
CUDA(1024ThreadsPerBlock):毫秒(每个样品经过的平均时间100样品)699.66.
以防万一,带32/ 的版本1024 nThreadsPerBlock是最快/最慢的版本.
我知道从主机到设备复制时有一种开销,反之亦然,但可能是因为我没有实现最快的代码.
因为我远不是一个CUDA专家:
我是否为此任务编写了最快的版本？我怎么能改进我的代码？我可以摆脱内核函数中的循环吗？
任何想法都赞赏.
虽然我描述了一个标准rowSum,但我对具有值的行的AND/ OR操作感兴趣(0;1},比如rowAND/ rowOR.也就是说,它不允许我利用cuBLAS乘法1的COL列向量技巧,正如一些评论家所建议的那样.
用户建议其他用户并在此认可:
忘记尝试编写自己的功能,使用Thrust库代替魔法来.
kan*_*yin 14
既然你提到你需要除了和之外的一般还原算法.我会尝试在这里提供3种方法.内核方法可能具有最高性能.推力方法最容易实现.cuBLAS方法仅适用于sum并且具有良好的性能.
这是一篇非常好的文档,介绍如何优化标准并行缩减.标准减少可分为两个阶段.
对于多减少(减少垫子的行)问题,只有阶段1就足够了.这个想法是每个线程块减少1行.有关每个线程块多行或每个多线程块1行的进一步考虑,可以参考@Novak提供的文章.这可以更好地改善性能,特别是对于形状不好的矩阵.
一般多次还原可以thrust::reduction_by_key在几分钟内完成.您可以在此处找到一些讨论.使用CUDA Thrust确定每个矩阵列中的最小元素及其位置.
但是thrust::reduction_by_key并不假设每行具有相同的长度,因此您将获得性能损失.另一篇文章如何在CUDA中规范化矩阵列并获得最大性能？给thrust::reduction_by_key出行和之间的分析比较和cuBLAS方法.它可以让您对性能有基本的了解.
矩阵A的行/列的总和可以看作矩阵向量乘法,其中向量的元素都是1.它可以用以下matlab代码表示.
y = A * ones(size(A,2),1);
哪里y是A的总和.
cuBLAS libary cublas<t>gemv()为此操作提供高性能矩阵向量乘法功能.
时序结果表明,该程序仅比简单读取A的所有元素慢10~50%,这可以看作是该操作性能的理论上限.
可以通过三种方式使用CUDA Thrust来减少矩阵的行数（它们可能不是唯一的，但解决这一点超出了范围）。同一个 OP 也认识到，对于此类问题，最好使用 CUDA Thrust。此外，使用cuBLAS的方法也是可能的。
方法#1 - reduce_by_key
这是Thrust 示例页面中建议的方法。它包括一个使用make_discard_iterator.
方法#2 - transform
这是CUDA Thrust的 Robert Crovella 建议的方法：reduce_by_key 仅针对数组中的某些值，基于“键”数组中的值。
方法#3 - inclusive_scan_by_key
这是 Eric 在How to normalize matrix columns in CUDA with max performance 中建议的方法？.
方法#4 - cublas<t>gemv
它用于cuBLAS gemv将相关矩阵乘以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，它比其他方法更方便，前提是创建计划所需的时间可以在计算过程中分摊。
| 归档时间： |
|
| 查看次数： |
10328 次 |
| 最近记录： |