RcppArmadillo очень удобно , и я использую его сам все время . Поскольку весь R cpp* код будет вызываться из R, мы можем предположить, что некоторые функциональные возможности присутствуют.
Это включает LAPACK и BLAS, и объясняет, почему мы можем использовать RcppArmadillo «без связывания», даже если в документации Armadillo ясно указано, что вам нужны LAPACK и BLAS. Почему? Хорошо , потому что R уже имеет LAPACK и BLAS .
(Кроме того, это приводит к значительным ранним проблемам, если и только если R был построен со своим собственным подмножеством LAPACK в виде некоторых сложных подпрограмм были недоступны. Баптист был сильно поражен этим, так как его пакеты нуждаются в этом (ред.). На протяжении многих лет Брайан Рипли был наиболее полезным в добавлении этих недостающих подпрограмм в LAPACK R. И никогда не возникало проблем, когда R создавался с внешний LAPACK и BLAS, как обычно для например, пакет Debian / Ubuntu, который я поддерживаю.)
Здесь вам нужен SuperLU. Это необязательно, и это ваша работа , чтобы убедиться, что это связано. Короче говоря, он не просто работает волшебным образом. И это трудно автоматизировать, так как это связано со связыванием, которое позволяет нам зависеть от платформы и требований к установке, которые мы не можем легко контролировать.
Но вопрос не нов, и на самом деле весь R cpp Галерея постов на топи c.
Редактировать: И с однострочным из этого поста, который подходит для моей системы, ваш код тоже работает отлично ( как только я исправлю 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)
*/