Rccp功能-SegFault错误和内存管理

Question

我需要遍历矩阵中的列，并对列索引较高的每一行中的所有值求和。

我已经使用for循环和rowSums很好地完成了工作，因为我熟悉基本的R代码。我写了一个C ++函数来加快运行时间-下面的代码。在RStudio中运行时，遇到致命错误。当我在R中运行时，由于未映射内存，因此会出现段错误。

当我运行例如10行的矩阵时不会发生错误-但理想情况下要运行10k。

我需要在某处分配内存吗？

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
Rcpp::NumericVector PVCalc_cpp(Rcpp::NumericMatrix x, int y) {  
int nrow = x.nrow(), ncol = x.ncol();
RNGScope scope;

Rcpp::NumericVector out(nrow);
double total = 0;

      for (int i = 0; i < nrow; i++) {
        total = 0;
        for (int j = y; j < ncol; j++){
          total += x(i, j+1);
        }
        out(i) =  floor((100. * total)+.5)/100;
      }
  return out;
}

Answer 1

问题是，当j = ncol - 1试图访问该行中的x(i, j+1)时>

total += x(i, j+1);
与尝试访问x(i, ncol)相同，并且您超出范围。如果我正确理解了您的问题，则只需将j+1正确传递即可将j更改为y。因此，我们可以将您的代码更改为以下内容：

#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::NumericVector PVCalc_cpp(const Rcpp::NumericMatrix& x, int y) {  
    int nrow = x.nrow(), ncol = x.ncol();
    Rcpp::NumericVector out(nrow);
    double total = 0;
    for (int i = 0; i < nrow; i++) {
        total = 0;
        for (int j = y; j < ncol; j++){
            total += x(i, j);
        }
        /* The following will always be equal to total but rounded;
           to type less R code for the R comparison,
           I omit the rounding here
         */
        // out(i) = floor((100. * total)+.5)/100;
        out(i) = total;
    }
    return out;
}
然后我们可以验证这是否可行，并且比rowSums()更快：

## Load required packages
library(Rcpp)
library(microbenchmark)
## Source C++ code
sourceCpp("so.cpp")
## Generate example data
x <- matrix(1:9, nrow = 3)
x
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
## Check results are the same
rowSums(x[ , -1])
[1] 11 13 15
PVCalc_cpp(x, 1)
[1] 11 13 15
## Benchmark small example
library(microbenchmark)
microbenchmark(base = rowSums(x[ , -1]),
               rcpp = PVCalc_cpp(x, 1))
Unit: microseconds
 expr   min     lq     mean median     uq      max neval
 base 6.285 6.7275  8.30454 7.0355 7.4375   54.616   100
 rcpp 2.648 2.8440 20.27161 3.3035 3.5520 1687.282   100
## Check larger example
set.seed(123)
x <- matrix(rnorm(1e6), nrow = 1e3)
y <- sample(seq_len(ncol(x)), size = 1)
all.equal(rowSums(x[ , -(1:y)]), PVCalc_cpp(x, y))
[1] TRUE
microbenchmark(base = rowSums(x[ , -(1:y)]),
               rcpp = PVCalc_cpp(x, y))
Unit: milliseconds
 expr      min       lq     mean   median       uq       max neval
 base 4.579241 5.179657 7.107276 6.588472 8.188014 14.179009   100
 rcpp 3.375648 3.654743 4.435746 4.454249 5.093247  6.024412   100