使用 CUDA 推力确定最小元素及其在每个矩阵列中的位置
Determining the least element and its position in each matrix column with CUDA Thrust
我有一个相当简单的问题，但我无法找到一个优雅的解决方案。
我有一个推力代码，它可以生成包含值的相同大小c向量。假设这些c向量中的每一个都有一个索引。我希望每个向量位置都获得值最低的c向量的索引：
例：
C0 = (0,10,20,3,40)
C1 = (1,2 ,3 ,5,10)
结果，我会得到一个包含具有最小值的C向量的索引的向量：
result = (0,1 ,1 ,0,1)
我曾考虑过使用 thrust zip 迭代器来做到这一点，但遇到了问题：我可以压缩所有c向量并实现一个任意转换，该转换接受元组并返回其最低值的索引，但是：
- 如何迭代元组的内容？
- 据我了解，元组最多只能存储
10个元素，而且可以有10c个向量。
然后我想过这样做：而不是c单独的向量，将它们全部附加到单个向量C中，然后生成引用位置的键并按键执行稳定的排序，这将从同一位置将向量条目重新组合在一起。在示例中，将给出：
C = (0,10,20,3,40,1,2,3,5,10)
keys = (0,1 ,2 ,3,4 ,0,1,2,3,4 )
after stable sort by key:
output = (0,1,10,2,20,3,3,5,40,10)
keys = (0,0,1 ,1,2 ,2,3,3,4 ,4 )
然后生成具有向量中位置的键，使用c向量的索引压缩输出，然后使用自定义函子执行 reduce by 键，对于每个归约输出具有最小值的索引。在示例中：
input = (0,1,10,2,20,3,3,5,40,10)
indexes= (0,1,0 ,1,0 ,1,0,1,0 ,1)
keys = (0,0,1 ,1,2 ,2,3,3,4 ,4)
after reduce by keys on zipped input and indexes:
output = (0,1,1,0,1)
但是，如何为按键操作减少编写这样的函子？
因为向量的长度必须相同。最好将它们连接在一起并将它们视为矩阵 C。
然后，您的问题就变成了在行主矩阵中查找每列的 min 元素的索引。可以按如下方式解决。
- 将行专业更改为副专业;
- 查找每列的索引。
在步骤 1 中，您建议使用 stable_sort_by_key 重新排列元素顺序，这不是一种有效的方法。由于在给定矩阵的 #row 和 #col 的情况下，可以直接计算重排。在推力中，可以使用排列迭代器完成，如下所示：
thrust::make_permutation_iterator(
c.begin(),
thrust::make_transform_iterator(
thrust::make_counting_iterator((int) 0),
(_1 % row) * col + _1 / row)
)
在第 2 步中，reduce_by_key可以完全按照您的要求进行操作。在您的情况下，归约二进制操作函子很容易，因为已经定义了元组(压缩向量的元素(上的比较来比较元组的第一个元素，并且推力支持它为
thrust::minimum< thrust::tuple<float, int> >()
整个程序如下所示。推力 1.6.0+ 是必需的，因为我在花哨的迭代器中使用占位符。
#include <iterator>
#include <algorithm>
#include <thrust/device_vector.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/iterator/discard_iterator.h>
#include <thrust/reduce.h>
#include <thrust/functional.h>
using namespace thrust::placeholders;
int main()
{
const int row = 2;
const int col = 5;
float initc[] =
{ 0, 10, 20, 3, 40, 1, 2, 3, 5, 10 };
thrust::device_vector<float> c(initc, initc + row * col);
thrust::device_vector<float> minval(col);
thrust::device_vector<int> minidx(col);
thrust::reduce_by_key(
thrust::make_transform_iterator(
thrust::make_counting_iterator((int) 0),
_1 / row),
thrust::make_transform_iterator(
thrust::make_counting_iterator((int) 0),
_1 / row) + row * col,
thrust::make_zip_iterator(
thrust::make_tuple(
thrust::make_permutation_iterator(
c.begin(),
thrust::make_transform_iterator(
thrust::make_counting_iterator((int) 0), (_1 % row) * col + _1 / row)),
thrust::make_transform_iterator(
thrust::make_counting_iterator((int) 0), _1 % row))),
thrust::make_discard_iterator(),
thrust::make_zip_iterator(
thrust::make_tuple(
minval.begin(),
minidx.begin())),
thrust::equal_to<int>(),
thrust::minimum<thrust::tuple<float, int> >()
);
std::copy(minidx.begin(), minidx.end(), std::ostream_iterator<int>(std::cout, " "));
std::cout << std::endl;
return 0;
}
剩下的两个问题可能会影响性能。
- 必须输出最小值，这不是必需的;
-
reduce_by_key是为具有可变长度的段设计的，但它可能不是对具有相同长度的段进行缩减的最快算法。
编写自己的内核可能是实现最高性能的最佳解决方案。
一个可能的想法，在这里建立矢量化排序的想法
-
假设我有这样的向量：
values: C = ( 0,10,20, 3,40, 1, 2, 3, 5,10) keys: K = ( 0, 1, 2, 3, 4, 0, 1, 2, 3, 4) segments: S = ( 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) -
将 K 和 S 压缩在一起以创建 KS
-
stable_sort_by_key使用 C 作为键，使用 KS 作为值：
stable_sort_by_key(C.begin(), C.end(), KS_begin); -
将重新排序的 C 和 K 向量压缩在一起，以创建 CK
-
stable_sort_by_key使用重新排序的 S 作为键，将 CK 作为值：
stable_sort_by_key(S.begin(), S.end(), CK_begin); -
使用排列迭代器或跨距范围迭代器访问新重新排序的 K 向量的每个第 N 个元素(0、N、2N、...(，以检索每个段中最小元素索引的向量，其中 N 是段的长度。
我还没有真正实现这个，现在只是一个想法。 也许由于我尚未观察到的某种原因，它不起作用。
segments(S(和keys(K(实际上是行索引和列索引。
对我来说似乎很奇怪，因为你的标题提到了"查找最大值的索引"，但你的大部分问题似乎都指的是"最低值"。 无论如何，通过更改我的算法的步骤 6，您可以找到任一值。
我很
好奇地测试以前的方法中哪一种更快。因此，我在下面的代码中实现了Robert Crovella的想法，为了完整起见，该代码还报告了Eric的方法。
#include <iterator>
#include <algorithm>
#include <thrust/random.h>
#include <thrust/device_vector.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/iterator/discard_iterator.h>
#include <thrust/reduce.h>
#include <thrust/functional.h>
#include <thrust/sort.h>
#include "TimingGPU.cuh"
using namespace thrust::placeholders;
template <typename Iterator>
class strided_range
{
public:
typedef typename thrust::iterator_difference<Iterator>::type difference_type;
struct stride_functor : public thrust::unary_function<difference_type,difference_type>
{
difference_type stride;
stride_functor(difference_type stride)
: stride(stride) {}
__host__ __device__
difference_type operator()(const difference_type& i) const
{
return stride * i;
}
};
typedef typename thrust::counting_iterator<difference_type> CountingIterator;
typedef typename thrust::transform_iterator<stride_functor, CountingIterator> TransformIterator;
typedef typename thrust::permutation_iterator<Iterator,TransformIterator> PermutationIterator;
// type of the strided_range iterator
typedef PermutationIterator iterator;
// construct strided_range for the range [first,last)
strided_range(Iterator first, Iterator last, difference_type stride)
: first(first), last(last), stride(stride) {}
iterator begin(void) const
{
return PermutationIterator(first, TransformIterator(CountingIterator(0), stride_functor(stride)));
}
iterator end(void) const
{
return begin() + ((last - first) + (stride - 1)) / stride;
}
protected:
Iterator first;
Iterator last;
difference_type stride;
};
/**************************************************************/
/* CONVERT LINEAR INDEX TO ROW INDEX - NEEDED FOR APPROACH #1 */
/**************************************************************/
template< typename T >
struct mod_functor {
__host__ __device__ T operator()(T a, T b) { return a % b; }
};
/********/
/* MAIN */
/********/
int main()
{
/***********************/
/* SETTING THE PROBLEM */
/***********************/
const int Nrows = 200;
const int Ncols = 200;
// --- 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 NR. 1 */
/******************/
timerGPU.StartCounter();
thrust::device_vector<float> d_min_values(Ncols);
thrust::device_vector<int> d_min_indices_1(Ncols);
thrust::reduce_by_key(
thrust::make_transform_iterator(
thrust::make_counting_iterator((int) 0),
_1 / Nrows),
thrust::make_transform_iterator(
thrust::make_counting_iterator((int) 0),
_1 / Nrows) + Nrows * Ncols,
thrust::make_zip_iterator(
thrust::make_tuple(
thrust::make_permutation_iterator(
d_matrix.begin(),
thrust::make_transform_iterator(
thrust::make_counting_iterator((int) 0), (_1 % Nrows) * Ncols + _1 / Nrows)),
thrust::make_transform_iterator(
thrust::make_counting_iterator((int) 0), _1 % Nrows))),
thrust::make_discard_iterator(),
thrust::make_zip_iterator(
thrust::make_tuple(
d_min_values.begin(),
d_min_indices_1.begin())),
thrust::equal_to<int>(),
thrust::minimum<thrust::tuple<float, int> >()
);
printf("Timing for approach #1 = %fn", timerGPU.GetCounter());
/******************/
/* APPROACH NR. 2 */
/******************/
timerGPU.StartCounter();
// --- Computing row indices vector
thrust::device_vector<int> d_row_indices(Nrows * Ncols);
thrust::transform(thrust::make_counting_iterator(0), thrust::make_counting_iterator(Nrows * Ncols), thrust::make_constant_iterator(Ncols), d_row_indices.begin(), thrust::divides<int>() );
// --- Computing column indices vector
thrust::device_vector<int> d_column_indices(Nrows * Ncols);
thrust::transform(thrust::make_counting_iterator(0), thrust::make_counting_iterator(Nrows * Ncols), thrust::make_constant_iterator(Ncols), d_column_indices.begin(), mod_functor<int>());
// --- int and float iterators
typedef thrust::device_vector<int>::iterator IntIterator;
typedef thrust::device_vector<float>::iterator FloatIterator;
// --- Relevant tuples of int and float iterators
typedef thrust::tuple<IntIterator, IntIterator> IteratorTuple1;
typedef thrust::tuple<FloatIterator, IntIterator> IteratorTuple2;
// --- zip_iterator of the relevant tuples
typedef thrust::zip_iterator<IteratorTuple1> ZipIterator1;
typedef thrust::zip_iterator<IteratorTuple2> ZipIterator2;
// --- zip_iterator creation
ZipIterator1 iter1(thrust::make_tuple(d_column_indices.begin(), d_row_indices.begin()));
thrust::stable_sort_by_key(d_matrix.begin(), d_matrix.end(), iter1);
ZipIterator2 iter2(thrust::make_tuple(d_matrix.begin(), d_row_indices.begin()));
thrust::stable_sort_by_key(d_column_indices.begin(), d_column_indices.end(), iter2);
typedef thrust::device_vector<int>::iterator Iterator;
// --- Strided access to the sorted array
strided_range<Iterator> d_min_indices_2(d_row_indices.begin(), d_row_indices.end(), Nrows);
printf("Timing for approach #2 = %fn", timerGPU.GetCounter());
printf("nn");
std::copy(d_min_indices_2.begin(), d_min_indices_2.end(), std::ostream_iterator<int>(std::cout, " "));
std::cout << std::endl;
return 0;
}
针对2000x2000大小的矩阵测试这两种方法，这是开普勒 K20c 卡上的结果：
Eric's : 8.4s
Robert Crovella's : 33.4s
相关文章：
- 将值指定给向量(2D）的向量中的某个位置
- OpenMP阵列性能较差
- 使用Unreal C++获取VR耳机的世界位置/方向
- 写入位置0x0000000C时发生访问冲突
- 如何将两个不同矢量的同一位置的两个元素组合在一起
- 访问特定阵列位置/索引时出现分段错误
- 访问字符阵列中不可用的内存位置(超出范围值）
- 如何将X重置回阵列中的原始位置
- mpi_bcast，将阵列缓冲区广播到接收缓冲区中的特定位置
- C 2D char阵列有时会在随机位置复制数据
- MPI-如何将杂志发送到阵列中的特定位置
- 查找阵列的位置-C
- 阵列相邻位置交换程序需要提示
- 检查阵列位置在C++中是否为空的 CPU 高效方法
- C++阵列在内存中的存储位置
- 从随机起始位置开始的螺旋阵列
- 检查阵列位置是否为空/空
- 检查阵列位置
- 阵列内存位置
- 检查二维井字阵列位置c++
