使用带有 RcppArmadillo 的 SuperLU 稀疏求解器答案

【问题标题】：Using SuperLU sparse solver with RcppArmadillo使用带有 RcppArmadillo 的 SuperLU 稀疏求解器
【发布时间】：2021-09-13 05:55:55
【问题描述】：

我正在尝试通过 RcppArmadillo 使用来自犰狳 (http://arma.sourceforge.net/docs.html#spsolve) 的 SparseLU 求解器：

#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。

【问题讨论】：

标签： r rcpp armadillo rcpparmadillo

【解决方案1】：

RcppArmadillo 超级方便，我自己一直在使用它。因为 all Rcpp* 代码将从 R 中调用，我们可以假设存在一些功能。

这包括 LAPACK 和 BLAS，并解释了为什么即使 Armadillo 文档明确声明您需要 LAPACK 和 BLAS，我们也可以“不链接”使用 RcppArmadillo。为什么？好吧因为 R 已经有 LAPACK 和 BLAS。

（顺便说一句，当且仅当 R 是用它自己的 LAPACK 子集构建时，这会导致相当多的早期问题，因为某些复杂的值例程不可用。Baptiste 受到了相当大的打击，因为他的包需要（ed ) 这些。多年来，Brian Ripley 在将那些缺失的例程添加到 R 的 LAPACK 中最有帮助。当使用外部 LAPACK 和 BLAS 构建 R 时，就像 eg Debian 一样，从来没有遇到过问题/Ubuntu 包我维护。）

这里你想要 SuperLU。这是可选的，你的工作是确保它被链接。简而言之，它确实不只是神奇地工作。而且很难实现自动化，因为它涉及到链接，这让我们无法轻松控制平台依赖和安装要求。

但这个问题并不新鲜，其实有an entire Rcpp Gallery post这个话题。

编辑： 使用该帖子中适合我的系统的单行代码，您的代码也可以正常工作（一旦我更正 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)
*/

【讨论】：