使用 CUDA Thrust 确定每个矩阵列中的最小元素及其位置
Posted
技术标签:
【中文标题】使用 CUDA Thrust 确定每个矩阵列中的最小元素及其位置【英文标题】：Determining the least element and its position in each matrix column with CUDA Thrust 【发布时间】：2013-07-15 22:59:10 【问题描述】：我有一个相当简单的问题，但我想不出一个优雅的解决方案。
我有一个推力代码，它产生包含值的相同大小的c 向量。假设每个c 向量都有一个索引。我想为每个向量位置获取值最低的c 向量的索引：
例子：
C0 = (0,10,20,3,40)
C1 = (1,2 ,3 ,5,10)
我会得到一个向量，其中包含具有最小值的C 向量的索引：
result = (0,1 ,1 ,0,1)
我曾考虑过使用推力 zip 迭代器，但遇到了问题：我可以压缩所有 c 向量并实现任意转换，该转换采用一个元组并返回其最小值的索引，但是：
-
如何遍历元组的内容？
据我了解，元组最多只能存储
10 元素，并且可以有更多的 10 c 向量。
然后我考虑过这样做：不要让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 向量的索引压缩输出，然后使用自定义函子按键执行归约，每次归约输出具有最低值的索引。在示例中：
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)
但是，如何为reduce by key操作编写这样的函子？
【问题讨论】：
您实际上是在尝试在行主矩阵中查找每列的最小元素的索引。 【参考方案1】：一个可能的想法，基于矢量化排序想法here
假设我有这样的向量：
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);
使用 permutation iterator 或 strided range iterator 访问新重新排序的 K 向量的每个第 N 个元素 (0, N, 2N, ...)，以检索每个段中的最小元素，其中 N 是段的长度。
我还没有真正实现这个，现在这只是一个想法。也许由于我尚未观察到的某种原因它不会起作用。
segments (S) 和 keys (K) 实际上是行和列索引。
您的问题对我来说似乎很奇怪，因为您的标题提到了“查找最大值索引”，但您的大部分问题似乎都指的是“最低值”。无论如何，通过更改我的算法的第 6 步，您可以找到任一值。
【讨论】：
感谢您的回复。它确实有效，除了在第 4 步和第 5 步中，C 和 S 应该被压缩并在 K 作为键上执行排序。你对标题是对的，我编辑了它:)【参考方案2】：因为向量的长度必须相同。最好将它们连接在一起，并将它们视为矩阵C。
那么您的问题就变成了在行主矩阵中查找每列的最小元素的索引。可以如下解决。
-
将 row-major 更改为 col-major；
查找每列的索引。
在步骤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> >()
整个程序如下所示。 Thrust 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 专为具有不同长度的段而设计，它可能不是减少相同长度段的最快算法。
编写自己的内核可能是获得最高性能的最佳解决方案。
【讨论】：
看来您应该能够使用另一个discard_iterator 来忽略minval 输出。
@JaredHoberock 我试过但无法使用 cuda5 + v1.6/v1.7 进行编译。一个错误？ error: no suitable conversion function from "thrust::detail::any_assign" to "float" exists【参考方案3】：
我很想测试以前哪种方法更快。因此，我在下面的代码中实现了 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 = %f\n", 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 = %f\n", timerGPU.GetCounter());
printf("\n\n");
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 大小的矩阵测试这两种方法，这是在 Kepler K20c 卡上的结果：
Eric's : 8.4s
Robert Crovella's : 33.4s
【讨论】：
以上是关于使用 CUDA Thrust 确定每个矩阵列中的最小元素及其位置的主要内容，如果未能解决你的问题，请参考以下文章