与 Rcpp 并行化下的索引错误

Question

假设我必须对 A 和 B 进行矩阵运算。我想用 OpenMP 编写一些 RcppArmadillo 代码来创建一个矩阵，其中 3 列和行等于 A 的列数乘以 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)
*/

我该如何解决这个问题？

Answer 1

要调试这样的问题，尽可能简化问题很重要。

在这种情况下，这意味着：

删除并行化。
减小函数的输入大小。
- 10 比 100 更容易看到。
为变量值添加跟踪语句。
运行代码

主要问题在于 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;
            ~~~~~                   ^~~~~

运行ning 代码然后给出：

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

与 Rcpp 并行化下的索引错误

Indexing error under parallelization with Rcpp

c++

r

openmp

rcpp