'Rcpp' integration for 'GNU GSL' vectors and matrices The 'GNU Scientific Library' (or 'GSL') is a collection of numerical routines for scientific computing. It is particularly useful for C and C++ programs as it provides a standard C interface to a wide range of mathematical routines. There are over 1000 functions in total with an extensive test suite. The 'RcppGSL' package provides an easy-to-use interface between 'GSL' data structures and R using concepts from 'Rcpp' which is itself a package that eases the interfaces between R and C++. This package also serves as a prime example of how to build a package that uses 'Rcpp' to connect to another third-party library. The 'autoconf' script, 'inline' plugin and example package can all be used as a stanza to write a similar package against another library.
This package uses Rcpp to connect the R system to the GNU GSL, a collection of numerical routines for scientific computing, particularly its vector and matrix classes.
lm() for OLS regressionThe fastLm() function included as file src/fastLm.cpp in the package:
#include <RcppGSL.h> #include <gsl/gsl_multifit.h>#include <cmath> // [[Rcpp::export]]Rcpp::List fastLm(const RcppGSL::Matrix &X, const RcppGSL::Vector &y) { int n = X.nrow(), k = X.ncol(); double chisq; RcppGSL::Vector coef(k); // to hold the coefficient vector RcppGSL::Matrix cov(k,k); // and the covariance matrix // the actual fit requires working memory we allocate and free gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (n, k); gsl_multifit_linear (X, y, coef, cov, &chisq, work); gsl_multifit_linear_free (work); // assign diagonal to a vector, then take square roots to get std.error Rcpp::NumericVector std_err; std_err = gsl_matrix_diagonal(cov); // need two step decl. and assignment std_err = Rcpp::sqrt(std_err); // sqrt() is an Rcpp sugar function return Rcpp::List::create(Rcpp::Named("coefficients") = coef, Rcpp::Named("stderr") = std_err, Rcpp::Named("df.residual") = n - k); } This example comes from the complete example package included in RcppGSL
and is from the file inst/examples/RcppGSLExample/src/colNorm.cpp
#include <RcppGSL.h>#include <gsl/gsl_matrix.h>#include <gsl/gsl_blas.h> // [[Rcpp::export]]Rcpp::NumericVector colNorm(const RcppGSL::Matrix & G) { int k = G.ncol(); Rcpp::NumericVector n(k); // to store results for (int j = 0; j < k; j++) { RcppGSL::VectorView colview = gsl_matrix_const_column (G, j); n[j] = gsl_blas_dnrm2(colview); } return n; // return vector} On CRAN, here and on its package page.
Dirk Eddelbuettel and Romain Francois
GPL (>= 2)