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.
The 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
}
- GNU GSL library (eg libgsl0-dev on Debian or Ubuntu)
- Rcpp for seamless R and C++ integration
On CRAN and here.
Dirk Eddelbuettel and Romain Francois
GPL (>= 2)