c - 通过 R 在 C 中操作矩阵

Question

我想用C写一些矩阵操作函数，然后传给R，R提供矩阵，得到操作结果。我有一个如下所示的测试函数(请不要介意它做了什么，在我原来的函数中，我需要从每一行中选择一个随机元素，并对它们进行一些计算，然后返回一个由这些随机组成的数组从每一行中选择元素，换句话说，我必须有两个 for 循环来遍历矩阵的所有元素)。

void multMat(double **A, int *r, int *c, double *s)
{
    int i, j;

    for (i = 0; i < *r; ++i)
    {

       for (j = 0; j < *c; ++j) 
       {
           if (j == 5)
               s[i] = A[i][j] * A[i][0];
       }
    }
}

我用 R CMD SHLIB multMat.c 编译了它，它为我生成了 multMat.so。然后，在 R 端，我有这样的东西:

dyn.load("multMat.so")

multMat <- function(A)
{
  .C("multMat", A=as.double(A), r=as.integer(nrow(A)), c=as.integer(ncol(A)), s=as.double(nrow(A)))
}

之后，我像在 R Studio 中一样创建了一个测试矩阵并调用了这个函数:

A <- matrix(1:100, 10, 10)
multMat(A)

问题是当我运行这个函数时，R Studio 崩溃了。我想 C 函数的定义方式存在一些问题。有什么想法吗？

Answer 1

哪里错了

您将 A 定义为 double **,

void multMat(double **A, int *r, int *c, double *s)

同时传递一个double *:

.C("multMat", A=as.double(A), r=as.integer(nrow(A)), c=as.integer(ncol(A)), s=as.double(nrow(A)))

您应该使用一维数组重写您的 C 函数。将您的函数定义为:

void multMat(double *A, int *r, int *c, double *s)

并将 A[i][j] 替换为 A[j * r + i]。 (如果我没有误会，r 是主要维度。)

性能问题:

目前是以j-loop为内层循环的i-j循环，所以你在最内层循环中扫描矩阵的一行。这缓存不友好。您应该交换循环以获得 j-i 循环。

我想您实际上已经意识到了缓存问题。在 C 中，矩阵按行优先顺序存储，因此 i-j 循环是最优的；但在 R 中，矩阵按列优先顺序存储，因此 j-i 循环是最优的。

也许矩阵的不同存储方式会导致您现有代码出现一些问题。你可能会在这里三思而后行。您的原始 C 代码假定它采用按行存储的矩阵，而如果您在 R 中初始化矩阵并将其提供给 C，则它按列存储。可能需要进行一些更改。如果这导致您的代码更改太多，您可以尝试将 R 矩阵的转置传递给 C。

您还应该在您的 C 代码中使用更多的本地/自动变量，而不是使用指针。例如替换

for (i = 0; i < *r; ++i)

通过

int r_local = *r;
for (i = 0; i < r_local; ++i)

您可以通过CPU 寄存器重用 和指令减少 获得性能提升(无需在每次迭代时取消引用该点)。

R 支持的 C/FORTRAN 数据类型

So, there is no way to directly pass the matrix from R, and then use the casual convention for matrix, meaning A[i][j] in C side?

不，没有。 R 不支持 double **。 R中的矩阵以column-major-format存储为long vector，属于double *类型。

R storage mode  C type        FORTRAN type
logical         int *         INTEGER
integer         int *         INTEGER
double          double *      DOUBLE PRECISION
complex         Rcomplex *    DOUBLE COMPLEX
character       char **       CHARACTER*255
raw unsigned    char *        none

参见 section 5.2 of Writing R extensions .