关于在R中使用C程序的一些问题

R是一个优秀的统计计算语言,但是因为它是解释型语言, 所以在对数组元素的迭代运算方面会很慢。 在R用C语言程序可以既保留R的易用性又可以在必要时提高速度。 本文讲述在Windows环境下如何用Borland的C编译器来完成R和C的结合。

假设我们要用C编码的问题是两个向量的卷积问题, 当然,R中已经有convolve可以实现, 我们这里只是作为一个例子。 两个无穷向量x和y的卷积定义为z[i] = sum_{k} x[i-k]*y[k], 两个有限长度向量的卷积是把向量两端填零后卷积然后取非填零部分计算得到的结果。 C程序为:

void convolve(double *a, int *na, double *b, int *nb, double *ab){
    int i, j, nab = *na + *nb - 1;
    for(i = 0; i < nab; i++)
    ab[i] = 0.0;
    for(i = 0; i < *na; i++)
    for(j = 0; j < *nb; j++)
    ab[i + j] += a[i] * b[j];
}
这里所有参数都使用指针, 是因为R中基本的运算单元是向量, 即C中的数组(指针)。

一、如何用Borland自由C++编译器生成Win32环境下的DLL(无图形界面成分)

首先要获得Borland的自由编译器,网址为 http://www.borland.com/bcppbuilder/freecompiler/ 假设安装在c:\borland\bcc55中。 在bin子目录中,生成一个bcc32.cfg文件, 包含如下行∶

-I"c:\Borland\Bcc55\include"
-L"c:\Borland\Bcc55\lib"
生成一个ilink32.cfg,包含如下行∶
-L"c:\Borland\Bcc55\lib"
生成如下的批处理文件mkdll.bat:
REM Make DLL for R using borland free compiler 5.5
PATH=c:\borland\bcc55\bin;C:\;c:\windows;c:\windows\command
bcc32 -O2 -WDE %1

假设我们有一个C程序文件testdll.c, 其中函数convolve是我们希望输出的, 我们需要另外生成一个文件testdll.def, 其中包含如下行∶

LIBRARY testdll
EXPORTS 
 _convolve

其中函数名前的下划线是编译器加的。 如果有多个输出函数可以写在后面,每个函数名占一行。

然后,在DOS窗口中运行

mkdll testdll.c
即可得到所需的testdll.dll库。 注意mkdll.bat文件应该在路径中, 否则应该使用全路径。

二、如何用BorlandCBuilder生成Win32环境下的DLL(无图形界面成分)

从文件菜单选择new...,然后选择Console wizard, 选Window type为console,Execution type为EXE, 即可生成程序框架, 把要输出到DLL中的函数写到源程序文件中, 并且要输出的函数用

extern  "C"  __declspec(dllexport) __stdcall
定义,如
extern  "C"  __declspec(dllexport) __stdcall 
void convolve(double *a, int *na, double *b, int *nb, double *ab){
。编译后可以得到所需的DLL库。

三、如何为R的Windows版准备用.C界面调用C或C++的动态链接库

用如上方法得到得到DLL库testdll.dll后, 进入R,运行

    dyn.load("testdll.dll")
可以把DLL库调入,注意如果库不在当前工作路径可以使用全路径, 如"c:/work/testdll.dll"。 然后,在R中定义如下包装函数conv:
conv <- function(a, b){
    .C("_convolve", as.double(a), as.integer(length(a)),
       as.double(b), as.integer(length(b)),
       ab=double(length(a)+length(b)-1))$ab
}

运行如:

    > conv(c(1,2), c(10,20,30))
    [1] 10 40 70 60

注意:

(1) R中所有对象传递到C中都是数组(即指针), 所以在我们的C程序中虽然参量na是一个整数标量我们还是用指针表示。

(2) R函数".C"的第一参数是C中的函数名, 用Borland 5.5 编译时会自动加上前导下划线, 其它参数都必须用as.double,as.int等包裹起来, 否则无法与C中的类型对应。".C"返回一个列表, 其中包含每一个调用参数(除函数名外), 有名调用的参数结果列表中元素名就使用给定的名字。 如本例中".C"函数返回一个包含元素"$ab"和四个无名元素的列表。

(3) 如果修改了C程序,要重新调入DLL库,可以先用dyn.unload("test.dll") 卸载然后再用dyn.load调入。

四、如何为R的Windows版准备用.Call界面调用C或C++的动态链接库

用.C界面可以很容易地从R中调用一个字符终端格式地Windows DLL, 但是不能实现从C中直接使用R对象的功能。 如果想从C中直接使用R对象, 就需要使用R的.Call界面来调用动态链接入的函数, 而在C程序中,使用SEXP来说明R对象, 并嵌入头文件R.h和Rdefines.h。 具体操作步骤如下:

(1) 安装Borland 自由C++编译器并设置好。

(2) 获得R的源程序包并安装在适当位置,如C:\R-1.2.1。

(3) 安装R Windows运行版本于适当位置,如C:\rw1021。

(4) 进入DOS命令窗口,确保程序搜索路径中有Borland编译器的路径, 进入C:\rw1021\bin,运行

     IMPDEF R.DEF R.DLL
这样可以得到R.DLL库中所有输出函数列表于R.DEF中。 然后,运行
     IMPLIB R.LIB R.DEF
这样可以得到用来与R.DLL链接的库R.LIB。

(5) 生成一个用来制作DLL库的批命令文件mkrdll.bat,内容为

REM Make DLL calling R using borland free compiler 5.5
PATH=c:\borland\bcc55\bin;C:\;c:\windows;c:\windows\command
bcc32 -u -O2 -WDE -IC:\R-1.2.1\src\include -LC:\rw1021\bin %1 R.lib

(6) 以上步骤只需要做一次。对C的源文件call.c,只要运行 mkrdll call.c 就可以生成所需得到DLL库call.dll。以卷积为例,C源程序如下:

#include 
#include 

SEXP convolve2(SEXP a, SEXP b)
{
    int i, j, na, nb, nab;
    double *xa, *xb, *xab;
    SEXP ab;
    PROTECT(a = AS_NUMERIC(a));
    PROTECT(b = AS_NUMERIC(b));
    na = LENGTH(a); nb = LENGTH(b); nab = na + nb - 1;
    PROTECT(ab = NEW_NUMERIC(nab));
    xa = NUMERIC_POINTER(a); xb = NUMERIC_POINTER(b);
    xab = NUMERIC_POINTER(ab);
    for(i = 0; i < nab; i++) xab[i] = 0.0;
    for(i = 0; i < na; i++)
    for(j = 0; j < nb; j++) xab[i + j] += xa[i] * xb[j];
    UNPROTECT(3);
    return(ab);
}

(7) 在R中,用

    > dyn.load("call.dll")
调入DLL库,然后用如
    > .Call("_convolve", c(1,2), c(10,20,30))
    [1] 10 40 70 60
用R内建函数计算为
    > convolve(c(1,2), rev(c(10,20,30)), type="open")
    [1] 10 40 70 60

注:

(1) 关于如何在C中调用R对象及函数的问题详见R随机手册 “Writing R Extensions”的第4章。

(2) 从上面的例子中我们可以看到, 为了能够从C中使用R对象,需要调入头函数R.h和Rdefines.h。 所有R对象都要说明成SEXP类型。

(3) 可以用AS_NUMERIC(a)把R对象a强制转换为数值向量类型的R对象, 以确保输入参数类型是正确的,但这一步最好是在R包裹代码中完成。 用NEW_NUMERIC(nab)可以生成一个新的R数值向量对象,长度nab。

(4) 所有在C代码中生成或重新赋值的R对象如上面的a,b,ab在赋值时都要用 PROTECT宏保护起来,格式如PROTECT(ab = NEW_NUMERIC(nab))。 在推出函数之前要用UNPROTEC函数去掉对新对象的保护, 参数为保护的对象的个数。

(5) 为了得到R数值向量的元素,用xa=NUMERIC_POINTER(a)可以得到该向量的头指针, 这是一个普通C的double类型指针,我们就可以把xa当作一个C数组处理了。

(6) 函数可以返回一个R对象,即返回值类型为SEXP。