使用RcppArmadillo与SuperLU稀疏求解器

Question

我尝试使用Armadillo的SparseLU求解器(http://arma.sourceforge.net/docs.html#spsolve)，通过RcppArmadillo实现:

#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)
*/

我遇到了错误undefined reference to 'superlu_free'，我猜测可能是缺少库链接。 你有什么思路解决这个问题吗？

---

我使用的是Windows 10操作系统。

- F. Privé

Answer 1

RcppArmadillo非常方便，我自己也一直在使用。因为所有的Rcpp*代码都将从R中调用，我们可以假设某些功能已经存在。这包括LAPACK和BLAS，这就解释了为什么我们即使没有链接Armadillo文档中明确指出需要LAPACK和BLAS也可以使用RcppArmadillo。

顺便提一下，如果R是使用其自己的LAPACK子集构建的，则会出现一些问题，因为某些复杂值例程不可用。多年来，Brian Ripley添加了这些缺失的例程到R的LAPACK中，这对Baptiste非常有帮助。如果R是使用外部LAPACK和BLAS构建的，则不会出现任何问题，这在我维护的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)
*/

- Dirk Eddelbuettel