[R-pkg-devel] How to obtain intercept of intercept-only glm in Fortran?
Ivan Krylov
kry|ov@r00t @end|ng |rom gm@||@com
Sat May 11 14:23:20 CEST 2019
On Fri, 10 May 2019 16:17:42 +0000
"Wang, Zhu" <wangz1 using uthscsa.edu> wrote:
> Are there any examples or links for me to follow through more closely?
Calling R functions from C++ is described at
<http://dirk.eddelbuettel.com/code/rcpp/Rcpp-quickref.pdf> and
elsewhere in Rcpp documentation. An example follows:
--------------8<--------------glmfit.cpp--------------8<--------------
#include <algorithm>
#include <Rcpp.h>
using namespace Rcpp;
extern "C" double intercept_glm(size_t n, const double * response) {
// access functions from default environment
Function glm_fit("glm.fit"), coef("coef");
// intercept-only model: response ~ 1
NumericVector x(n);
x.fill(1);
// I couldn't find a way to wrap a double* into a NumericVector
// without copying anything, sorry; perhaps someone else
// can offer a solution
NumericVector y(n);
std::copy_n(response, n, y.begin());
// call the R function, convert the result back
return as<double>(coef(glm_fit(x, y)));
}
--------------8<--------------glmfit.cpp--------------8<--------------
Since this function is extern "C" and uses only primitive C types, it
should be fairly easy to call from Fortran. (C is the lingua franca of
programming languages). Fortran-C interoperability is well described in
"Modern Fortran Explained" by Metcalf et al. Here is the Fortran side
of the code:
--------------8<--------------callglm.f90--------------8<--------------
subroutine callglm(ret)
use, intrinsic :: iso_c_binding, only: c_size_t, c_double
! using iso_c_binding here
! - to get correct type of ret when R calls the function
! - to convert variables before calling C function
implicit none
! using F77-style arguments to match expectations of .Fortran()
real(c_double), intent(out) :: ret
! toy data to compare against R code later
real :: y(10) = [10, 11, 20, 9, 10, 8, 11, 45, 2, 3]
! the interface block declares an extern "C" function
interface
! double intercept_glm(size_t n, const double * response)
function intercept_glm(n, response) bind(c)
use, intrinsic :: iso_c_binding
real(c_double) :: intercept_glm
integer(c_size_t), value :: n
real(c_double) :: response(*)
end function
end interface
! call the function as you would call any other function
ret = intercept_glm(int(size(y), c_size_t), real(y, c_double))
end subroutine
--------------8<--------------callglm.f90--------------8<--------------
For a quick test, make sure that you have Rcpp installed and run:
# adjust R version and path if your library is elsewhere
PKG_CPPFLAGS='-g -I ~/R/x86_64-pc-linux-gnu-library/3.3/Rcpp/include' \
R CMD SHLIB callglm.f90 glmfit.cpp
R
library(Rcpp)
dyn.load('callglm.so') # change extension if needed
.Fortran('callglm', ret=numeric(1))
# $ret
# [1] 12.9
coef(glm.fit(rep(1, 10), c(10, 11, 20, 9, 10, 8, 11, 45, 2, 3)))
# [1] 12.9
To use this in a package, place both files in the src/ subdirectory of
your package and add LinkingTo: Rcpp in the DESCRIPTION.
--
Best regards,
Ivan
More information about the R-package-devel
mailing list