# Using the GSL to compute eigenvalues

[This article was first published on

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

**Rcpp Gallery**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Two posts showed how to compute eigenvalues using Armadillo and using Eigen. As we also looked at using the

GNU GSL, this post will show how to conpute eigenvalues using GSL.

As mentioned in the previous GSL post, we instantiate C language pointers suitable for GSL (here the matrix `M`

). Those *must* be freed manually, as shown before the `return`

statement.

// [[Rcpp::depends(RcppGSL)]] #include <RcppGSL.h> #include <gsl/gsl_matrix.h> #include <gsl/gsl_eigen.h> // [[Rcpp::export]] Rcpp::NumericVector getEigenValues(Rcpp::NumericMatrix sM) { RcppGSL::matrix<double> M(sM); // create gsl data structures from SEXP int k = M.ncol(); RcppGSL::vector<double> ev(k); // instead of gsl_vector_alloc(k); gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc(k); gsl_eigen_symm (M, ev, w); gsl_eigen_symm_free (w); // we need to invoke wrap() here to help the compiler, but at least we save a loop // we also need to create a new Rcpp object as the RcppGSL one needs to be free'd. Rcpp::NumericVector res(Rcpp::wrap(ev)); M.free(); // important: GSL wrappers use C structure ev.free(); return res; // return results vector }

We can illustrate this easily via a random sample matrix.

set.seed(42) X <- matrix(rnorm(4*4), 4, 4) Z <- X %*% t(X) getEigenValues(Z) [1] 14.2100 2.4099 1.6856 0.3319

In comparison, R gets the same results (in reverse order) and also returns the eigenvectors.

eigen(Z) $values [1] 14.2100 2.4099 1.6856 0.3319 $vectors [,1] [,2] [,3] [,4] [1,] 0.69988 -0.55799 0.4458 -0.00627 [2,] -0.06833 -0.08433 0.0157 0.99397 [3,] 0.44100 -0.15334 -0.8838 0.03127 [4,] 0.55769 0.81118 0.1413 0.10493

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 about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.