Ошибка индексации при распараллеливании с Rcpp - PullRequest
Новая
       19

Ошибка индексации при распараллеливании с Rcpp

0 голосов
/ 16 февраля 2019

Предположим, мне нужно матрицу A и B. Я хочу написать некоторый код RcppArmadillo с OpenMP, который создает матрицу с 3 столбцами и строками, равными количеству столбцов A, умноженному на количество строк B.

Я написал этот код, но он вылетает, когда я пытаюсь его запустить.Я думаю, что я делаю ошибку при создании переменной row, но я не уверен, как ее исправить. 

#include "RcppArmadillo.h"
#include <omp.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;

// [[Rcpp::plugins(openmp)]]
// [[Rcpp::export]]
arma::mat my_matrix(const arma::mat & A, 
                    const arma::mat B) {

  const int nObservations = A.n_cols;
  const int nDraws = B.n_rows;
  const int nRows = nObservations*nRows;
  arma::mat out(nRows,3);
  int i,n,iter,row;
  omp_set_num_threads(2);
  #pragma omp parallel for private(i, n, iter, row)
  for(i = 0; i < nDraws; i++){
    for(n = 0; n < nObservations; n++) {
      row = i * nObservations + n ;
      out(row,0) = i+1 ;
      out(row,1) = n+1 ;
      out(row,2) = row+1 ;
    }
  }

  return out;
}

/*** R
set.seed(9782)
A <- matrix(rnorm(100), ncol = 5)
B <- matrix(rnorm(100), nrow = 10)


test <- my_matrix(A = A, B = B)
*/

Как я могу это исправить?

1 Ответ

0 голосов
/ 16 февраля 2019

Для устранения подобных проблем важно максимально упростить проблему.

В этом случае это означает:

  1. Удалить распараллеливание.
  2. Уменьшите входной размер до функции.
    • 10 намного легче увидеть, чем 100.
  3. Добавить операторы трассировки для значений переменных.
  4. Код выполнения

Основная проблема заключается в том, как строится nRows: 

const int nRows = nObservations * nRows;
                               // ^^^^^ Self-reference

Переключите его на: 

const int nRows = nObservations * nDraws;

Затем добавьте параллелизацию, и все будет хорошо.

Пример упрощенного кода с операторами трассировки для отладки. 

#include "RcppArmadillo.h"
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;

// [[Rcpp::export]]
arma::mat my_matrix(const arma::mat & A, 
                    const arma::mat B) {

  const int nObservations = A.n_cols;
  const int nDraws = B.n_rows;
  const int nRows = nObservations * nRows;

  // Show initialization information
  Rcpp::Rcout << "nObservations: " << nObservations << std::endl 
              << "nDraws: " << nDraws << std::endl 
              << "nRows: " << nRows << std::endl;

  arma::mat out(nRows, 3);

  // Show trace of matrix construction
  Rcpp::Rcout << "out - rows: " << out.n_rows << std::endl 
              << "out - columns: " << out.n_cols << std::endl;

  int i, n, iter, row;
  for(i = 0; i < nDraws; ++i){
    for(n = 0; n < nObservations; ++n) {
      row = i * nObservations + n;
      // Show trace statement of index being accessed
      Rcpp::Rcout << "Output row access id: " << row << std::endl;

      out(row, 0) = i + 1;
      out(row, 1) = n + 1;
      out(row, 2) = row + 1;
    }
  }

  return out;
}

Компиляция этого фрагмента кода дает два предупреждения, относящиеся к неиспользуемым переменным ... 

file69cab2726a1.cpp:13:13: warning: unused variable 'iter' [-Wunused-variable]
  int i, n, iter, row;
            ^
file69cab2726a1.cpp:11:37: warning: variable 'nRows' is uninitialized when used within its own initialization [-Wuninitialized]
  const int nRows = nObservations * nRows;
            ~~~~~                   ^~~~~

Запуск кода затем дает: 

set.seed(9782)
A <- matrix(rnorm(10), ncol = 5)
B <- matrix(rnorm(10), nrow = 10)

test <- my_matrix(A = A, B = B)
# nObservations: 5
# nDraws: 10
# nRows: 0
# out - rows: 0
# out - columns: 3
# Output row access id: 0
# 
# error: Mat::operator(): index out of bounds
...