[R] Use R function in C code

Douglas Bates bates at stat.wisc.edu
Thu Apr 15 15:07:52 CEST 2004

James Wettenhall <wettenhall at wehi.edu.au> writes:

> Hi,
> On Fri Apr 2, xt_wang wrote:
> > I want to use R function Matrix inverse in my c code, please 
> > tell me how I can.
> > If there is a sample which can tell me how it works. It will 
> > be fantastic.
> A good place to start learning how to interface R with C is the 
> "Writing R Extensions" manual installed locally, or:
> http://cran.r-project.org/doc/manuals/R-exts.pdf
> http://rweb.stat.umn.edu/R/doc/manual/R-exts.html
> See also the article "In Search of C/C++ & Fortran Routines" in: 
> http://cran.r-project.org/doc/Rnews/Rnews_2001-3.pdf
> I guess you'd have to decide whether you want to compile R as a 
> shared library or not.  I've never done this.
> But do you really want to call R from C to get a matrix inverse?  
> Can't you just use a CLAPACK routine?
> http://www.netlib.org/clapack/
> I'm not sure which CLAPACK subroutine is best for your purposes, 
> but I have in the past used dgels (double-precision gaussian 
> elimination and least-squares) and found it to be good.
> I assume you have asked the question: 
> "Do I really need an inverse?"

Exactly.  The fact that we write a formula using an inverse of a
matrix does not mean that this is a good way to compute the result.

Versions 0.8-2 and higher of the Matrix package contain a vignette
called Comparisons.pdf that provides comparative timings of several
different approaches to solving a linear least squares problem.

> If Ax = b,
> x = inv(A) * b 
> is computationally about twice as expensive as Gaussian 
> Elimination.

Actually the "ge" in "dgels" stands for "general", not Gaussian
elimination.  For square coefficient matrices the 'solve' function in
R now uses an LU decomposition (dgetrf) followed by dgetrs for
non-missing "b" or dgetri for the inverse ("b" missing).

The general rule is that if you are solving a specific set of linear
equations use dgetrf and dgetrs.  If you decide that you really do
need an inverse, which is surprisingly rare, use dgetrf and dgetri.

More information about the R-help mailing list