@@ -86,17 +86,18 @@ sourceCpp(code='
8686// [[Rcpp::export]]
8787Rcpp::List fun(Rcpp::NumericVector Yr, Rcpp::NumericMatrix Xr){
8888
89- int i,j, n = Xr.nrow(), k = Xr.ncol();
89+ int i, j, n = Xr.nrow(), k = Xr.ncol();
9090 double chisq;
9191
92- gsl_matrix *X = gsl_matrix_alloc (n, k);
93- gsl_vector *y = gsl_vector_alloc (n);
94- gsl_vector *c = gsl_vector_alloc (k);
95- gsl_matrix *cov = gsl_matrix_alloc (k, k);
92+ RcppGSL::Matrix X(n, k); // allocate a gsl_matrix<double> of dim n, k
93+ RcppGSL::Vector y(n); // allocate a gsl_vector<double> of length n
94+ RcppGSL::Vector c(k); // allocate a gsl_vector<double> of length k
95+ RcppGSL::Matrix cov(k, k); // allocate a gsl_matrix<double> of dim k, k
96+
9697 for (i = 0; i < n; i++) {
9798 for (j = 0; j < k; j++)
98- gsl_matrix_set (X, i, j, Xr(i,j) );
99- gsl_vector_set (y, i, Yr(i));
99+ X( i, j) = Xr(i, j );
100+ y[i] = Yr(i); // Note vector requires [] not ()
100101 }
101102
102103 gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (n, k);
@@ -105,17 +106,13 @@ Rcpp::List fun(Rcpp::NumericVector Yr, Rcpp::NumericMatrix Xr){
105106
106107 Rcpp::NumericVector coefr(k), stderrestr(k);
107108 for (i = 0; i < k; i++) {
108- coefr(i) = gsl_vector_get(c,i) ;
109- stderrestr(i) = sqrt(gsl_matrix_get( cov, i,i));
109+ coefr(i) = c[i] ;
110+ stderrestr(i) = sqrt(cov( i,i));
110111 }
111- gsl_matrix_free (X);
112- gsl_vector_free (y);
113- gsl_vector_free (c);
114- gsl_matrix_free (cov);
115112
116113
117- return Rcpp::List::create( Rcpp::Named( "coef", coefr) ,
118- Rcpp::Named( "stderr", stderrestr) );
114+ return Rcpp::List::create( Rcpp::Named("coef") = coefr,
115+ Rcpp::Named("stderr") = stderrestr);
119116}' )
120117fun
121118}
0 commit comments