Indexing error under parallelization with Rcpp

Question

Suppose I have to matrix A and B. I want to write some RcppArmadillo code with OpenMP that creates a matrix with 3 columns and rows equal to the number of columns of A times the number of rows of B.

I wrote this code but it crashes when I try to run it. My best guess is that I'm making an error when creating the row variable but I'm not sure how to fix it.

#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)
*/

How can I fix this?

Answer 1

To debug problems like this, it's important to simplify the problem as much as possible.

In this case, that means:

1. Remove parallelization.
2. Lower the input size to the function.
   • 10 is much easier to see than 100.
3. Add trace statements for variable values.
4. Run code

The main issue is with how nRows is being constructed:

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

Switch it to:

const int nRows = nObservations * nDraws;

Then re-add the parallelization and all should be well.

---

Example of simplified code with trace statements for debugging.

#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;
}

Compiling this piece of code gives two warnings related to unused variables...

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;
            ~~~~~                   ^~~~~

Running the code then gives:

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