RcppArmadillo超级方便,我自己一直都在用。因为所有Rcpp* 代码都将从 R 中调用,所以我们可以假设存在一些功能。
这包括 LAPACK 和 BLAS,并解释了为什么即使 Armadillo 文档明确声明您需要 LAPACK 和 BLAS,我们也可以“不链接”使用 RcppArmadillo。为什么?好吧,因为 R 已经有 LAPACK 和 BLAS。
(顺便说一句,当且仅当 R 是用它自己的 LAPACK 子集构建时,这会导致相当多的早期问题,因为某些复杂的值例程不可用。Baptiste 受到了相当大的打击,因为他的包需要(ed)这些。多年来,Brian Ripley 在将那些缺失的例程添加到 R 的 LAPACK 中最有帮助。当使用外部 LAPACK 和 BLAS 构建 R 时,从来没有遇到过问题,例如我维护的 Debian/Ubuntu 软件包。)
在这里,您需要 SuperLU。这是可选的,您的工作是确保将其链接起来。简而言之,它不仅仅是神奇地工作。而且很难自动化,因为它涉及链接,这使我们无法轻松控制平台依赖和安装要求。
但这个问题并不新鲜,事实上有一个关于这个主题的整个 Rcpp Gallery 帖子。
编辑:并且使用适合我的系统的该帖子的单行代码,您的代码也可以正常工作(一旦我更正了Rcpp::depends
它需要的双括号:
R> Sys.setenv("PKG_LIBS"="-lsuperlu")
R> sourceCpp("answer.cpp")
R> library(Matrix)
R> K <- sparseMatrix(i = c(1, 2, 1), j = c(1, 2, 2), x = c(1, 1, 0.5), symmetric = TRUE)
R> x <- runif(2)
R> sp_solver(K, x)
[,1]
[1,] 0.264929
[2,] 0.546050
R>
更正的代码
#define ARMA_USE_SUPERLU
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
arma::vec sp_solver(arma::sp_mat K, arma::vec x) {
arma::superlu_opts opts;
opts.symmetric = true;
arma::vec res;
arma::spsolve(res, K, x, "superlu", opts);
return res;
}
/*** R
library(Matrix)
K <- sparseMatrix(i = c(1, 2, 1), j = c(1, 2, 2), x = c(1, 1, 0.5), symmetric = TRUE)
x <- runif(2)
sp_solver(K, x)
*/