rcpp/html/lmGSL_8R_source.html

00001 #
00002 # lm() via GNU GSL -- building on something from the Intro to HPC tutorials
00003 #
00004 # Copyright (C) 2010 Dirk Eddelbuettel and Romain Francois
00005 #
00006 # This file is part of Rcpp.
00007 #
00008 # Rcpp is free software: you can redistribute it and/or modify it
00009 # under the terms of the GNU General Public License as published by
00010 # the Free Software Foundation, either version 2 of the License, or
00011 # (at your option) any later version.
00012 #
00013 # Rcpp is distributed in the hope that it will be useful, but
00014 # WITHOUT ANY WARRANTY; without even the implied warranty of
00015 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00016 # GNU General Public License for more details.
00017 #
00018 # You should have received a copy of the GNU General Public License
00019 # along with Rcpp.  If not, see <http://www.gnu.org/licenses/>.
00020
00021 suppressMessages(require(Rcpp))
00022 suppressMessages(require(inline))
00023
00024 lmGSL <- function() {
00025
00026     src <- '
00027
00028     Rcpp::NumericVector Yr(Ysexp);
00029     Rcpp::NumericMatrix Xr(Xsexp);
00030
00031     int i,j,n = Xr.nrow(), k = Xr.ncol();
00032     double chisq;
00033
00034     gsl_matrix *X = gsl_matrix_alloc (n, k);
00035     gsl_vector *y = gsl_vector_alloc (n);
00036     gsl_vector *c = gsl_vector_alloc (k);
00037     gsl_matrix *cov = gsl_matrix_alloc (k, k);
00038     for (i = 0; i < n; i++) {
00039         for (j = 0; j < k; j++)
00040             gsl_matrix_set (X, i, j, Xr(i,j));
00041         gsl_vector_set (y, i, Yr(i));
00042     }
00043
00044     gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (n, k);
00045     gsl_multifit_linear (X, y, c, cov, &chisq, work);
00046     gsl_multifit_linear_free (work);
00047
00048     Rcpp::NumericVector coefr(k), stderrestr(k);
00049     for (i = 0; i < k; i++) {
00050         coefr(i) = gsl_vector_get(c,i);
00051         stderrestr(i) = sqrt(gsl_matrix_get(cov,i,i));
00052     }
00053     gsl_matrix_free (X);
00054     gsl_vector_free (y);
00055     gsl_vector_free (c);
00056     gsl_matrix_free (cov);
00057
00058
00059     return Rcpp::List::create( Rcpp::Named( "coef", coefr),
00060                                Rcpp::Named( "stderr", stderrestr));
00061     '
00062
00063     ## turn into a function that R can call
00064     ## compileargs redundant on Debian/Ubuntu as gsl headers are found anyway
00065     fun <- cfunction(signature(Ysexp="numeric", Xsexp="numeric"),
00066                      src,
00067                      includes="#include <gsl/gsl_multifit.h>",
00068                      Rcpp=TRUE,
00069                      cppargs="-I/usr/include",
00070                      libargs="-lgsl -lgslcblas")
00071 }