与Rcpp并行化下的索引错误

Question

假设我必须使用矩阵A和B.我想用OpenMP编写一些RcppArmadillo代码，它创建一个矩阵，其中3列和行的行数等于B的行数的A倍。

我写了这段代码，但是当我尝试运行它时崩溃了。我最好的猜测是我在创建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

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

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

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