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
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
