[R] Use R function in C code

[Douglas Bates]

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.