一个朴素的 C 解决方案是
library(inline)
fun4 <-
cfunction(c(x="numeric", a="numeric", b="numeric"), body4,
language="C")
body4 <- "
R_len_t len = Rf_length(x);
SEXP result = Rf_allocVector(REALSXP, len);
const double aa = REAL(a)[0], bb = REAL(b)[0], *xp = REAL(x);
double *rp = REAL(result);
for (int i = 0; i < len; ++i)
if (xp[i] < aa)
rp[i] = aa;
else if (xp[i] > bb)
rp[i] = bb;
else
rp[i] = xp[i];
return result;
"
fun4 <-
cfunction(c(x="numeric", a="numeric", b="numeric"), body4,
language="C")
使用简单的并行版本(正如 Dirk 指出的那样,这是CFLAGS = -fopenmp
在 ~/.R/Makevars 中,并且在支持 openmp 的平台/编译器上)
body5 <- "
R_len_t len = Rf_length(x);
const double aa = REAL(a)[0], bb = REAL(b)[0], *xp = REAL(x);
SEXP result = Rf_allocVector(REALSXP, len);
double *rp = REAL(result);
#pragma omp parallel for
for (int i = 0; i < len; ++i)
if (xp[i] < aa)
rp[i] = aa;
else if (xp[i] > bb)
rp[i] = bb;
else
rp[i] = xp[i];
return result;
"
fun5 <-
cfunction(c(x="numeric", a="numeric", b="numeric"), body5,
language="C")
和基准
> z <- runif(1e7)
> benchmark(fun1(z,0.25,0.75), fun4(z, .25, .75), fun5(z, .25, .75),
+ replications=10)
test replications elapsed relative user.self sys.self
1 fun1(z, 0.25, 0.75) 10 9.087 14.609325 8.335 0.739
2 fun4(z, 0.25, 0.75) 10 1.505 2.419614 1.305 0.198
3 fun5(z, 0.25, 0.75) 10 0.622 1.000000 2.156 0.320
user.child sys.child
1 0 0
2 0 0
3 0 0
> identical(res1 <- fun1(z,0.25,0.75), fun4(z,0.25,0.75))
[1] TRUE
> identical(res1, fun5(z, 0.25, 0.75))
[1] TRUE
在我的四核笔记本电脑上。假设数字输入、无错误检查、NA 处理等。