Rcpp函数 - SegFault错误& 内存管理

Question

我需要对矩阵中的列进行循环，并对列指数较高的每一行中的所有值进行求和。

我已经用for循环和rowSums做得很好，没有问题，因为我熟悉基本的R代码。我写了一个C++函数来加快运行速度--代码如下。当我在RStudio中运行时，我遇到了一个致命的错误。当我在R中运行时，我得到一个segfault错误，因为内存没有被映射。

当我运行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.
           Notice also the use of [] access operator instead of ();
           this is more efficient.
         */
        // 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 5.591 5.9210  8.61073  6.475 6.786  137.125   100
 rcpp 2.337 2.5795 19.90118  3.035 3.222 1651.094   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 5.377342 6.052347 6.954338 6.482907 7.834190 11.580706   100
 rcpp 1.447596 1.909504 2.085185 2.023343 2.158256  3.159366   100