Можно ли (когда-либо) использовать R cpp :: IntegerVector с OpenMP? - PullRequest
Новая
       101

Можно ли (когда-либо) использовать R cpp :: IntegerVector с OpenMP?

2 голосов
/ 07 мая 2020

Я, вероятно, жаден до производительности, но я наблюдал значительный прирост производительности при объединении R cpp и OpenMP, возможно, незаконным образом. Я понимаю, что «Вызов любого R API из многопоточного кода является« только для экспертов »и категорически не рекомендуется». но я не совсем понимаю, когда код C ++ может неявно делать это на R cpp векторов. Я понимаю, что RcppParallel имеет класс RVector, но я понимаю, что это может включать в себя снятие копии вектора, что может sh отменить любой прирост производительности.

Я заметил, что в галерее R cpp есть (https://gallery.rcpp.org/articles/gerber-statistic/), который, похоже, обращается к NumericMatrix HIST_RETURN_RAW внутри параллельного l oop, поэтому кажется, что «некоторый» доступ разрешен, но я знаю / верю, что некоторые оболочки (например, Rcpp::List) будут не допускаются. Является ли атомарность отличительной характеристикой c?

Конкретно, можно ли использовать OpenMP в следующих случаях? (т.е. все ли они потокобезопасны, совместимы с управлением памятью R и свободны от неопределенного поведения?) 

#include <Rcpp.h>
#ifdef _OPENMP
#include <omp.h>
#endif
using namespace Rcpp;


// 1. No 'outside' R vector, but IntegerVector created outside omp region

// [[Rcpp::export]]
IntegerVector fastInitalize(int n, int nThread = 1) {
  IntegerVector out = no_init(n);
#pragma omp parallel for num_threads(nThread)
  for (int i = 0; i < n; ++i) {
    out[i] = 0;
  }
  return out;
}


// 2. Simple access

// [[Rcpp::export]]
IntegerVector AddOMP(IntegerVector x, IntegerVector y, int nThread = 1) {
  R_xlen_t N = x.length();
  if (N != y.length()) {
    stop("Lengths differ");
  }
  IntegerVector out = no_init(N);

#pragma omp parallel for num_threads(nThread)
  for (R_xlen_t i = 0; i < N; ++i) {
    out[i] = x[i] + y[i];
  }
  return out;
}


// 3. Access, copy inside Rcpp

// [[Rcpp::export]]
IntegerVector pmax0xy(IntegerVector x, IntegerVector y, int nThread = 1) {
  R_xlen_t N = x.length();
  if (N != y.length()) {
    stop("Lengths differ");
  }
  IntegerVector out = clone(y);

#pragma omp parallel for num_threads(nThread)
  for (R_xlen_t i = 0; i < N; ++i) {
    if (x[i] > 0) {
      out[i] = 0;
    }
  }
  return out;
}


// 4. Access with omp + reduction

// [[Rcpp::export]]
int firstNonzero(IntegerVector x, int nThread = 1) {
  R_xlen_t N = x.length();
  int out = N;
#pragma omp parallel for num_threads(nThread) reduction(min : out)
  for (R_xlen_t i = 0; i < N; ++i) {
    if (x[i] != 0) {
      out = (i < out) ? i : out;
    }
  }
  return out;
}


// 5. Access with omp array reduction

// [[Rcpp::export]]
IntegerVector count_one_to_ten(IntegerVector x, int nThread = 1) {
  R_xlen_t N = x.length();
  if (N >= INT_MAX) {
    stop("Possibly too many numbers.");
  }
  const int TEN = 10;
  int numbers[TEN] = {}; // what if 10 was large?
#if defined _OPENMP && _OPENMP >= 201511
#pragma omp parallel for num_threads(nThread) reduction(+:numbers[:TEN])
#endif
  for (R_xlen_t i = 0; i < N; ++i) {
    int xi = x[i];
    if (xi < 1 || xi > 10) {
      continue;
    }
    numbers[xi - 1] += 1;
  }
  IntegerVector out(TEN);
  for (int n = 0; n < TEN; ++n) {
    out[n] = numbers[n];
  }
  return out;
}



// You can include R code blocks in C++ files processed with sourceCpp
// (useful for testing and development). The R code will be automatically
// run after the compilation.
//

/*** R
x <- sample(1:1200, size = 1e6, replace = TRUE)
y <- sample(1:1200, size = 1e6, replace = TRUE)

fastInitalize(1e6)[1]

head(AddOMP(x, y))
head(AddOMP(x, y, 2))

head(pmax0xy(x, y))
head(pmax0xy(x, y, 2))

firstNonzero(x)
firstNonzero(x, 2)

count_one_to_ten(x, 2)
*/
...