CUDA/Thrust: Как суммировать столбцы чередующегося массива?

Question

Используя Thrust, вы можете суммировать строки чередующегося массива (т. Е. Подкрепленного вектором), как показано в примере here.

То, что я хотел бы сделать, - это сумма столбцов .

Я попытался с помощью аналогичную конструкцию, т.е .:

// convert a linear index to a column index 
template <typename T> 
struct linear_index_to_col_index : public thrust::unary_function<T,T> 
{ 
    T C; // number of columns 

    __host__ __device__ 
    linear_index_to_col_index(T C) : C(C) {} 

    __host__ __device__ 
    T operator()(T i) 
    { 
    return i % C; 
    } 
}; 

// allocate storage for column sums and indices 
thrust::device_vector<int> col_sums(C); 
thrust::device_vector<int> col_indices(C); 

// 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_col_index<int>(C)), 
    thrust::make_transform_iterator(thrust::counting_iterator<int>(0), linear_index_to_col_index<int>(C)) + (R*C), 
    array.begin(), 
    col_indices.begin(), 
    col_sums.begin(), 
    thrust::equal_to<int>(), 
    thrust::plus<int>()); 

Однако это приводит только суммирования первого столбца, остальные игнорируются. Моя догадаться, почему это происходит, является то, что, как было отмечено в reduce_by_keydocs:

Для каждой группы последовательных ключей в диапазоне [keys_first, keys_last), которые равны, reduce_by_key копирует первый элемент из группы к keys_output. [Упор мое]

Если я правильно понимаю, потому что ключи в строке итератора являются последовательными (т.е. индексы [0 - (С-1)] дает 0, то [C - (2C- 1)] даст 1 и т. Д.), Они в конечном итоге суммируются.

Но итератор столбца будет отображать индексы [0 - (C-1)] на [0 - (C-1)], а затем снова начать, индексы [C - (2C-1)] будут отображаться на [ 0 - (C-1)] и т. Д., Делая значения, полученные не последовательными.

Для меня это неинтуитивно, я ожидал, что все точки данных, назначенные одному и тому же ключу, будут сгруппированы вместе, но это еще одно обсуждение.

В любом случае, мой вопрос: как можно суммировать столбцы чередующегося массива с помощью Thrust?

Answer 1

Эти операции (суммирование строк, суммирующих столбцов и т. Д.), Как правило, связаны с пропускной способностью на GPU. Поэтому мы можем рассмотреть, как построить алгоритм, который будет оптимально использовать пропускную способность памяти GPU. В частности, мы хотели бы, чтобы наши базовые обращения к памяти, созданные из кода тяги, были объединены, если это возможно. В двух словах это означает, что соседние потоки графического процессора будут считываться из соседних местоположений в памяти.

Оригинал row-summing example отображает это свойство: смежные потоки, порожденные тягой, будут считывать смежные элементы в памяти. Например, если у нас есть строки R, то мы увидим, что первые потоки R, созданные путем нажатия, будут считывать первую «строку» матрицы во время операции reduce_by_key. Поскольку ячейки памяти, связанные с первой строкой, сгруппированы вместе, мы получаем совместный доступ.

Одним из подходов к решению этой проблемы (как суммировать столбцы) было бы использование аналогичной стратегии для примера суммирования строк, но используйте permutation_iterator, чтобы заставить потоки, являющиеся частью одной и той же последовательности клавиш, читать a столбец данных вместо строка данных. Этот итератор перестановки будет принимать базовый массив, а также последовательность отображения.Эта последовательность сопоставлений создается с помощью transform_iterator с использованием special functor, примененного к counting_iterator, для преобразования индекса линейной (большой строки) в индекс столбца, так что первые потоки C будут считывать элементы первого столбца вместо первой строки. Поскольку первые потоки C будут принадлежать одной и той же последовательности клавиш, они будут суммироваться в операции reduce_by_key. Это то, что я называю методом 1 в коде ниже.

Однако этот метод имеет недостатки, связанные с тем, что смежные потоки больше не читают смежные значения в памяти - мы сломали коалесценцию, и, как мы увидим, влияние производительности заметно.

Для больших матриц, хранящихся в строке-мажорная порядке в памяти (упорядочение мы обсуждали в этой задаче), довольно оптимальный метод суммирования столбцов должен иметь каждый поток сумму отдельный столбец с для- петля. Это довольно просто реализовать в CUDA C, и мы можем аналогичным образом выполнить эту операцию в Thrust с соответствующим образом определенным функтором.

Я имею в виду это как метод 2 в приведенном ниже коде. Этот метод будет запускать столько потоков, сколько есть столбцов в матрице. Для матрицы с достаточно большим количеством столбцов (скажем, 10 000 или более) этот метод будет насыщать графический процессор и эффективно использовать доступную пропускную способность памяти. Если вы проверите функтора, вы увидите, что это несколько «необычная» адаптация тяги, но совершенно законная.

Вот код сравнения оба метода:

$ cat t994.cu 
#include <thrust/device_vector.h> 
#include <thrust/reduce.h> 
#include <thrust/iterator/permutation_iterator.h> 
#include <thrust/iterator/transform_iterator.h> 
#include <thrust/iterator/counting_iterator.h> 
#include <thrust/functional.h> 
#include <thrust/sequence.h> 
#include <thrust/transform.h> 

#include <iostream> 

#define NUMR 1000 
#define NUMC 20000 
#define TEST_VAL 1 

#include <time.h> 
#include <sys/time.h> 
#define USECPSEC 1000000ULL 

long long dtime_usec(unsigned long long start){ 

    timeval tv; 
    gettimeofday(&tv, 0); 
    return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start; 
} 


typedef int mytype; 

// from a linear (row-major) index, return column-major index 
struct rm2cm_idx_functor : public thrust::unary_function<int, int> 
{ 
    int r; 
    int c; 

    rm2cm_idx_functor(int _r, int _c) : r(_r), c(_c) {}; 

    __host__ __device__ 
    int operator() (int idx) { 
    unsigned my_r = idx/c; 
    unsigned my_c = idx%c; 
    return (my_c * r) + my_r; 
    } 
}; 


// convert a linear index to a column index 
template <typename T> 
struct linear_index_to_col_index : public thrust::unary_function<T,T> 
{ 
    T R; // number of rows 

    __host__ __device__ 
    linear_index_to_col_index(T R) : R(R) {} 

    __host__ __device__ 
    T operator()(T i) 
    { 
    return i/R; 
    } 
}; 

struct sum_functor 
{ 
    int R; 
    int C; 
    mytype *arr; 

    sum_functor(int _R, int _C, mytype *_arr) : R(_R), C(_C), arr(_arr) {}; 

    __host__ __device__ 
    mytype operator()(int myC){ 
    mytype sum = 0; 
     for (int i = 0; i < R; i++) sum += arr[i*C+myC]; 
    return sum; 
    } 
}; 



int main(){ 
    int C = NUMC; 
    int R = NUMR; 
    thrust::device_vector<mytype> array(R*C, TEST_VAL); 

// method 1: permutation iterator 

// allocate storage for column sums and indices 
    thrust::device_vector<mytype> col_sums(C); 
    thrust::device_vector<int> col_indices(C); 

// compute column sums by summing values with equal column indices 
    unsigned long long m1t = dtime_usec(0); 
    thrust::reduce_by_key(thrust::make_transform_iterator(thrust::counting_iterator<int>(0), linear_index_to_col_index<int>(R)), 
    thrust::make_transform_iterator(thrust::counting_iterator<int>(R*C), linear_index_to_col_index<int>(R)), 
    thrust::make_permutation_iterator(array.begin(), thrust::make_transform_iterator(thrust::make_counting_iterator<int>(0), rm2cm_idx_functor(R, C))), 
    col_indices.begin(), 
    col_sums.begin(), 
    thrust::equal_to<int>(), 
    thrust::plus<int>()); 
    cudaDeviceSynchronize(); 
    m1t = dtime_usec(m1t); 
    for (int i = 0; i < C; i++) 
    if (col_sums[i] != R*TEST_VAL) {std::cout << "method 1 mismatch at: " << i << " was: " << col_sums[i] << " should be: " << R*TEST_VAL << std::endl; return 1;} 
    std::cout << "Method1 time: " << m1t/(float)USECPSEC << "s" << std::endl; 

// method 2: column-summing functor 

    thrust::device_vector<mytype> fcol_sums(C); 
    thrust::sequence(fcol_sums.begin(), fcol_sums.end()); // start with column index 
    unsigned long long m2t = dtime_usec(0); 
    thrust::transform(fcol_sums.begin(), fcol_sums.end(), fcol_sums.begin(), sum_functor(R, C, thrust::raw_pointer_cast(array.data()))); 
    cudaDeviceSynchronize(); 
    m2t = dtime_usec(m2t); 
    for (int i = 0; i < C; i++) 
    if (fcol_sums[i] != R*TEST_VAL) {std::cout << "method 2 mismatch at: " << i << " was: " << fcol_sums[i] << " should be: " << R*TEST_VAL << std::endl; return 1;} 
    std::cout << "Method2 time: " << m2t/(float)USECPSEC << "s" << std::endl; 
    return 0; 
} 
$ nvcc -O3 -o t994 t994.cu 
$ ./t994 
Method1 time: 0.034817s 
Method2 time: 0.00082s 
$ 

Очевидно, что при достаточно большой матрицы, метод 2 существенно быстрее, чем метод 1.

Если вы не знакомы с подстановок итераторы, принимают посмотрите на thrust quick start guide.