Код C с openmp, вызванным из R, дает противоречивые результаты

Question

Ниже приведен фрагмент кода C, запускаемый из R, который используется для сравнения каждой строки матрицы с вектором.Число идентичных значений хранится в первом столбце матрицы из двух столбцов.
Я знаю, что это легко сделать в R (как это делается для проверки результатов), но это первый шаг для более сложного использования.case.
Когда openmp не используется, он работает нормально.Когда используется openmp, он дает коррелированные (0.99), но противоречивые результаты.
Вопрос1: Что я делаю не так?
Вопрос2: Я использую двойной цикл for длязаполните выходную матрицу (ret) нулями.Что было бы лучшим решением?

Кроме того, наблюдались несоответствия при использовании кода в пакете.Я попытался сделать код воспроизводимым, используя inline, но он не распознает операторы openmp (я попытался включить 'omp.h' в параметры cfunction, ...).
Вопрос3: Как мы можем заставить этот код работать с inline?

Я (тоже?) Далеко за пределами моей зоны комфорта по этой теме.

library(inline)

compare <- cfunction(c(x = "integer", vec = "integer"), "
  const int I = nrows(x), J = ncols(x);
  SEXP ret;
  PROTECT(ret = allocMatrix(INTSXP, I, 2));
  int *ptx = INTEGER(x), *ptvec = INTEGER(vec), *ptret = INTEGER(ret);

  for (int i=0; i<I; i++)
      for (int j=0; j<2; j++)
          ptret[j * I + i] = 0;

  int i, j;

  #pragma omp parallel for default(none) shared(ptx, ptvec, ptret) private(i,j)

  for (j=0; j<J; j++)
   for (i=0; i<I; i++)
      if (ptx[i + I * j] == ptvec[j]) {++ptret[i];}

  UNPROTECT(1);
  return ret;
")


N = 3e3
M = 1e4
m = matrix(sample(c(-1:1), N*M, replace = TRUE), nc = M)
v = sample(-1:1, M, replace = TRUE)

cc = compare(m, v)
cr = rowSums(t(t(m) == v))
all.equal(cc[,1], cr)

Answer 1

Благодаря комментариям выше я пересмотрел проблему гонки данных.
IIUC, мой цикл был распараллелен на j (столбцы).Затем каждый поток имел свое собственное значение i (строки), но возможные идентичные значения для потоков, которые затем пытались увеличить ptret[i] одновременно.Чтобы избежать этого, теперь я сначала включаю цикл i, так что только один поток будет увеличивать каждую строку.
Затем я понял, что могу переместить инициализацию нуля ptret в первом цикле.
Кажется, работает.Я получаю идентичные результаты, повышенную нагрузку на процессор и ускорение в 3-4 раза на своем ноутбуке.
Полагаю, это решает вопросы 1 и 2. Я более подробно рассмотрю проблему inline / openmp.
Код ниже, fwiw.

#include <omp.h>
#include <R.h>
#include <Rinternals.h>
#include <stdio.h>

SEXP c_compare(SEXP x, SEXP vec)
{
  const int I = nrows(x), J = ncols(x);
  SEXP ret;
  PROTECT(ret = allocMatrix(INTSXP, I, 2));
  int *ptx = INTEGER(x), *ptvec = INTEGER(vec), *ptret = INTEGER(ret);

  int i, j;
#pragma omp parallel for default(none) shared(ptx, ptvec, ptret) private(i, j)
  for (i = 0; i < I; i++) {
    // init ptret to zero
    ptret[i] = 0;
    ptret[I + i] = 0;

    for (j = 0; j < J; j++)
      if (ptx[i + I * j] == ptvec[j]) {
        ++ptret[i];
      }
    }

  UNPROTECT(1);
  return ret;
}