RcppParallel алгоритм быстрого векторного квадрата - PullRequest
Новая
       67

RcppParallel алгоритм быстрого векторного квадрата

0 голосов
/ 03 мая 2020

После руководства Я написал функцию, позволяющую подключить NumericVector через пакет RcppParallel. Код выглядит следующим образом: 

#include <RcppArmadillo.h>
#include <RcppParallel.h>

using namespace Rcpp;
using namespace RcppArmadillo;
using namespace RcppParallel;

double sqr(double x)
{
  return(x * x);
}

// Parallel pow of vectors
struct ParallelVectorPowStruct : public Worker
{
  // source matrix
  const RVector<double> input;

  // powers matrix
  const RVector<double> input_powers;

  // destination matrix
  RVector<double> output;

  // type of power (0 - general, 1 - square, 2 - square root)
  int pow_type;

  // initialize with source and destination
  ParallelVectorPowStruct(const NumericVector input, NumericVector output, NumericVector input_powers, int pow_type) 
    : input(input), output(output), input_powers(input_powers), pow_type(pow_type) {}

  // take the powers of the range of elements requested
  void operator()(std::size_t begin, std::size_t end) {

    if (pow_type == 0)
    {
      // perform pow operation
      std::transform(input.begin() + begin, 
                     input.begin() + end, 
                     input_powers.begin(),
                     output.begin() + begin, 
                     ::pow);
    }

    if (pow_type == 1)
    {
      std::transform(input.begin() + begin, 
                     input.begin() + end,
                     output.begin() + begin, 
                     sqr);
    }

    if (pow_type == 2)
    {
      // perform pow operation
      std::transform(input.begin() + begin, 
                     input.begin() + end,
                     output.begin() + begin, 
                     ::sqrt);
    }

  }

};

// [[Rcpp::export]]
NumericVector ParallelVectorPow(NumericVector x, double value) {

  int pow_type = 0;

  if (value == 0.5)
  {
    pow_type = 2;
  }

  if (value == 2)
  {
    pow_type = 1;
  }

  // allocate the output matrix
  NumericVector output(x.size());

  //  allocate the powers matrix
  NumericVector input_powers(x.size());
  std::fill(input_powers.begin(), input_powers.end(), value);

  // ParallelVectorPowStruct functor
  ParallelVectorPowStruct parallelVectorPowStruct(x, output, input_powers, pow_type);

  // call parallelFor to do the work
  parallelFor(0, x.length(), parallelVectorPowStruct);

  // return the output matrix
  return output;
}

Обратите внимание, что я добавил два дополнительных варианта вычисления (когда pow_type = 1 и когда pow_type = 2) для степени 0,5 (квадрат root) и для Степень 2 (квадрат). Моя мотивация для создания этого материала была следующей:

Сначала я проанализировал производительность кода без введения этих двух дополнительных типов возможностей и обнаружил, что распараллеленная версия работает намного быстрее (обычно примерно в 8 раз на моих 8- частота ядра 3800 МГц)

Например, следующий код: 

m <- runif(1000000, 1, 100)
library(rbenchmark)
res <- benchmark(m ^ 3,
                 ParallelVectorPow(m, 3),
                 order = "relative")
res[,1:4]

дает: 

                     test replications elapsed relative
2 ParallelVectorPow(m, 3)          100    0.45    1.000
1                     m^3          100    4.22    9.378

Однако все обстоит иначе, когда дело доходит до 0,5 и 2. Параллельная версия кода, используемая для работы намного медленнее, чем не распараллеленная для этих полномочий. Я решил проблему для мощности 0,5 с помощью функции :: sqrt, следуя подходу, упомянутому в руководстве, указанном выше (см. Код для pow_type == 2). Однако я не могу понять, как улучшить код для степени 2.

Так что рассмотрим следующий код: 

m <- runif(1000000, 1, 100)
library(rbenchmark)
res <- benchmark(m ^ 2,
                 ParallelVectorPow(m, 2),
                 order = "relative")
res[,1:4]

Это дает: 

                     test replications elapsed relative
1                     m^2          100    0.14    1.000
2 ParallelVectorPow(m, 2)          100    0.36    2.571

Можно ли ускорить его через фреймворк RcppParallel?

Будет очень здорово помочь!

PS

  1. Я думаю, что m ^ 2 и m ^ 0.5 делают специальные вызовы для некоторых быстрых c -языковых функций, в то время как m ^ other_power делает какую-то общую альтернативную подпрограмму.

  2. Я также проверил эту реализацию на различных подходах к квадрату на основе RcppArmadillo и обнаружил, что они медленнее, чем оба перечисленных выше .

...