C-Function 的 R-Wrapper 无法正常工作答案

【问题标题】：R-Wrapper for C-Function does not work properlyC-Function 的 R-Wrapper 无法正常工作
【发布时间】：2020-06-30 14:43:44
【问题描述】：

目前我正在向 C 介绍自己以扩展我的 R 函数。我用 C 语言编写了一个函数来进行一些计算，它们运行良好。但是，一旦我在 R 本身中编写了一个包装器，似乎就会出现一些错误。将我的 C 函数视为“colV”，将“abc”视为某个任意矩阵。

语句 (R) .Call("colV", abc, ncol(abc), nrow(abc)) 工作正常（每次，无论我使用它的频率如何），而

colV = function(x){
  nc = ncol(x)
  nr = nrow(x)
  .Call("colV", x, nc, nr)
}

在第三次使用时提供错误的结果：

> colV(abc)
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> colV(abc)
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> colV(abc)
[1] -8.370087e+22  6.254796e-01  6.774042e-01  1.709462e+00  7.250386e-01

如果我再次声明包装器，前两次运行正常，第三次尝试再次出现相同的结果。请注意，除了外观之外，实际上只有第一个值发生了变化！

有人知道包装器有什么问题吗？如前所述，如果我只使用 .Call 语句，则始终返回正确的结果：

> .Call("colV", abc, ncol(abc), nrow(abc))
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> .Call("colV", abc, ncol(abc), nrow(abc))
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> .Call("colV", abc, ncol(abc), nrow(abc))
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> .Call("colV", abc, ncol(abc), nrow(abc))
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> .Call("colV", abc, ncol(abc), nrow(abc))
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386
> .Call("colV", abc, ncol(abc), nrow(abc))
[1] 1.4274933 0.6254796 0.6774042 1.7094617 0.7250386

另外，我在使用 .C 接口时没有遇到过这个问题，但是由于速度和内存的原因，这对我来说不是一个解决方案。提前致谢艾琳

编辑：这是 C 代码：

#include <R.h>
#include <Rinternals.h>
#include <Rmath.h>


SEXP colV(SEXP y, SEXP n, SEXP r){
    int *nc = INTEGER(n);
    double *x = REAL(y);
    int d = length(y);
    int *nr = INTEGER(r);
    int i, j, z;
    //int d = nr * nc;
    double colMean[(*nc)];
    double xSq[(d)];
    double colMsq[(*nc)];
    double xSm[(*nc)];
    SEXP result;
    PROTECT(result = allocVector(REALSXP, (*nc)));
    memset(REAL(result), 0, (*nc) * sizeof(double));
    double *colVar = REAL(result);
    //PROTECT(colVar = NEW_DOUBLE(nc));
    int fr = ((*nr) - 1);


    for(z = 0; z < (d); z++){
        xSq[z] = x[z] * x[z];
    }

    for(i = 0; i < (*nc); i++){
        for(j = 0; j < (*nr); j++){
            colMean[i] += (x[(j + ((*nr) * i)) ]);
            xSm[i] += (xSq[(j + (*nr * i))]);
        }
        colMean[i] = (colMean[i] / (*nr));
        colMsq[i] = (*nr) * (colMean[i] * colMean[i]);
        colVar[i] = ((xSm[i] - colMsq[i]) / fr);
    }
    UNPROTECT(1);
    return(result);
}

【问题讨论】：

是否也可以看到 C 代码，因为没有它，这不是 MWE，因为它无法复制。
我可以提供它，我只是认为将其省略以提高可读性，因为直接调用 C 代码可以正常工作。但我会编辑它，没问题。
但是，对于不同的 c 代码，我遇到了完全相同的问题。只需调用 c 函数即可提供准确/预期的结果，而包装它会产生废话。所以，如果有人对这个问题有答案，请告诉我们原因。

标签： r c function call wrapper

【解决方案1】：

编辑：正如评论中所说，这（下）有效，但是对于另一个不同的 c-Function 会出现同样的问题。使用 .Call 直接提供预期的结果，而将 .Call 放在包装器中会出现问题。我还没找到原因。

我现在找到了问题的解决方案，但我不明白为什么它是解决方案。如果有人可以解释这一点，请随时这样做！

改变循环中 colMean 部分的计算改变了一切：

...
for(i = 0; i < (*nc); i++){
        for(j = 0; j < (*nr); j++){
            colMean[i] += ((x[(j + ((*nr) * i)) ]) / (*nr));
            xSm[i] += (xSq[(j + (*nr * i))]);
        }
        //colMean[i] = (colMean[i] / (*nr));
        colMsq[i] = (*nr) * (colMean[i] * colMean[i]);
        colVar[i] = ((xSm[i] - colMsq[i]) / fr);
    }
...

因此，将 nr 的除法直接移动到加法语句中即可解决问题。

祝大家有个美好的一天/晚上！

【讨论】：

到目前为止，我怀疑它与变量的范围有关。