[Rd] R + C + Lapack toy regression example

[Vinh Nguyen]

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" )
             , 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