Using Armadillo with SuperLU

February 17, 2017
By

(This article was first published on Rcpp Gallery, and kindly contributed to R-bloggers)

Armadillo is very versatile C++ library for linear algebra, brough
to R via the RcppArmadillo package. It
has proven to be very useful and popular, and is (as of February 2017) used by well over 300 CRAN
packages as indicated by the reverse depends / linking-to on its
CRAN page. Several earlier posts here on the
Rcpp Gallery also demonstrate different usage patterns.

Armadillo has a core focus on dense matrices, but continues to
expand its capabilities for sparse matrices. Basic operation are supported directly via the
templated header files, along with calls back into the default (dense) LAPACK and BLAS libraries we
can access easily as R uses them.

Armadillo also supports dedicated sparse solvers via the
SuperLU package. However, this requires access to
the SuperLU library. Many Linux distributions ship
it, see e.g. the Debian package page and the
Ubuntu package page; there is also a
homebrew recipe for OS X /
macOS (or other systems using brew). As of this writing, the version in the current Ubuntu
release is behind the version Debian. But it is the 5.2.* version that is in Debian that is also
required with current Armadillo versions 7.700.* so we prepared ‘backports’ via
Dirk’s PPA repo.)

Recently, a GitHub issue ticket asked how to
use SuperLU along with
RcppArmadillo. This post essentially
repeats the main answer, which was spread out over multiple posts, in a single article.

In a nutshell, two things need to happen:

  1. One needs to define the required variable ARMA_USE_SUPERLU which has to be done before the
    Armadillo headers are included. One possibility (shown below) is a #define statement right in the
    code file.

  2. One needs to tell the linker to use the SuperLU
    library. This step is of course not perfectly portable, and does require that the library be
    installed. (A genuinely portable solution would either need to test for presence of SuperLU, or include its
    sources. Both aspects are beyond the scope of this post._

Sys.setenv("PKG_LIBS"="-lsuperlu")

Now that R knows about this library, we can present the code requiring it:

// Important: this definition ensures Armadillo enables SuperLU
#define ARMA_USE_SUPERLU 1

#include 
// [[Rcpp::depends(RcppArmadillo)]]

using namespace arma;

// [[Rcpp::export]]
void superLuDemo() {
    sp_mat A = sprandu<sp_mat>(1000, 1000, 0.1);

    vec b = randu<vec>(1000);
    mat B = randu<mat>(1000, 5);

    vec x = spsolve(A, b);  // solve one system
    mat X = spsolve(A, B);  // solve several systems

    bool status = spsolve(x, A, b);  // use default solver
    if (status == false)  { Rcpp::Rcout << "no solution" << endl; }

    spsolve(x, A, b, "lapack");   // use LAPACK  solver
    spsolve(x, A, b, "superlu");  // use SuperLU solver

    Rcpp::Rcout << "Done.\n";
}

This code snippet defines a function superLuDemo() which we can call from R:

superLuDemo()
Done.

As the data generated here is random, we did not bother printing the (dense) result vector of
length 1000.

To leave a comment for the author, please follow the link and comment on their blog: Rcpp Gallery.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)