[Rd] R + C + Lapack toy regression example

[Douglas Bates]

On Wed, Sep 23, 2009 at 2:39 AM, Vinh Nguyen <vqnguyen at uci.edu> wrote:
> dear list,
>
> since matrix manipulations is often of interest in statistical
> computations, i'd like to get a working example of using Lapack for
> regression.  However, i run into an error.
>
> My matrix-lapack-example.c file:
> #include <R_ext/Lapack.h>
>
> void reg(const char* trans, const int* m, const int* n,
>         const int* nrhs, double* a, const int* lda,
>         double* b, const int* ldb,
>         double* work, const int* lwork, int* info)
> {
>  F77_CALL(dgels)(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info);
> }
>
> My matrix-lapack-example.R file:
> dyn.load( "matrix-lapack-example.so" )
>
> regression <- function( X, Y ){
>  m <- nrow( X )
>  n <- ncol( X )
>  stopifnot( m >= n , is.vector( Y ), n != length( Y ) )
>  mn <- min( m, n )
>  lwork <- mn + max( mn, 1 )
>  ### to be used with dgels
>  ### http://www.netlib.org/lapack/double/dgels.f
>  rslt <- .C( "reg"
>             , trans=as.character( "N" )

In your C code you are treating that object as a C character string
(char*) but that is not what is being passed.  Look at section 5.2 of
"Writing R Extensions" - a character variable in R is passed by .C as
a char**, not a char*.

I find it much easier to use the .Call interface instead of the .C
interface for applications like this.  You do need to learn some of
the mysteries of the functions declared in Rinternals.h but, once you
do, it is much easier to pass a matrix from R with and extract all the
fussy details within the C code.  Several of the C source code files
in the Matrix package are devoted to exactly the kind of operation you
want to carry out.  Look at the function lsq_dense_QR in
Matrix/src/dense.c, for example. (Although now that I look at it
myself I see some questionable programming practices - I should have
PROTECTed the result of coerceVector but it happens that it would not
have needed protection.  Rather than coercing I should just check
isInteger on the dim attribute.)
>             , m=as.integer( m ), n=as.integer( n )
>             , nrhs=as.integer( 1 ), a=as.double( X )
>             , lda=as.integer( m ), b=as.double( Y )
>             , ldb=as.integer( m ) , work=double( lwork )
>             , lwork=as.integer( lwork )
>             , info=integer( 1 ) )
>  ##F77_NAME(dgels)(const char* trans, const int* m, const int* n,
>  ##                const int* nrhs, double* a, const int* lda,
>  ##                double* b, const int* ldb,
>  ##                double* work, const int* lwork, int* info);
>  return( rslt )
> }
>
> n <- 100
> x <- rnorm( 100 )
> y <- 2.5 + 3*x + rnorm( n )
> X <- cbind( 1, x )
>
> temp <- regression( X=X, Y=y )
>
>
> dgels does linear least squares.  the C code compile fines with a
> warning (ld: warning, duplicate dylib
> /usr/local/lib/libgcc_s.1.dylib).  in R, i get the following when i
> run regression():
>> temp <- regression( X=X, Y=y )
> Parameter 1 to routine DGELS  was incorrect
> Mac OS BLAS parameter error in DGELS , parameter #0, (unavailable), is 0
>
> Process R exited abnormally with code 1 at Wed Sep 23 00:21:59 2009
>
> am i doing something wrong?  is my as.character() used wrong here?
>
> i know R uses fortran code for lm.  but suppose i wanted to do
> regression in C/C++.  is this lapack method the recommended / "best
> practice" route?  i just want to know whats the recommended method for
> doing matrix manipulations in C.  thanks!
>
> vinh
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>