On 20-04-2012, at 21:18, Jonathan Greenberg wrote:
> Ok, I figured out a solution and I'd like to get some feedback on this from
> the R-helpers as to how I could modify the following to be "package
> friendly" -- the main thing I'm worried about is how to dynamically set the
> "dyn.load" statement below correctly (obviously, its hard coded to my
> particular install, and would only work with windows since I'm using a
> .dll):
>
> Rdggev <- function(JOBVL=F,JOBVR=T,A,B)
> {
> # R implementation of the DGGEV LAPACK function (with generalized
> eigenvalue computation)
> # See http://www.netlib.org/lapack/double/dggev.f
> # coded by Jonathan A. Greenberg <j...@illinois.edu>
> dyn.load("C:\\Program Files\\R\\R-2.14.2\\bin\\x64\\Rlapack.dll")
>
> if(JOBVL)
> {
> JOBVL="V"
> } else
> {
> JOBVL="N"
> }
> if(JOBVR)
> {
> JOBVR="V"
> } else
> {
> JOBVR="N"
> }
> # Note, no error checking is performed.
> # Input parameters
> N=dim(A)[[1]]
> LDA=N
> LDB=N
> LDVL=N
> LDVR=N
> LWORK=as.integer(max(1,8*N))
> Rdggev_out <- .Fortran("dggev", JOBVL, JOBVR, N, A, LDA, B, LDB,
> double(N), double(N), double(N),
> array(data=0,dim=c(LDVL,N)), LDVL, array(data=0,dim=c(LDVR,N)), LDVR,
> double(max(1,LWORK)), LWORK, integer(1))
>
> names(Rdggev_out)=c("JOBVL","JOBVR","N","A","LDA","B","LDB","ALPHAR","ALPHAI",
> "BETA","VL","LDVL","VR","LDVR","WORK","LWORK","INFO")
>
> Rdggev_out$GENEIGENVALUES=(Rdggev_out$ALPHAR+Rdggev_out$ALPHAI)/Rdggev_out$BETA
> return(Rdggev_out)
> }
>
See this:
<start R code>
# This works on Mac OS X
# Change as needed for other systems
# or compile geigen into a standalone shared object.
dyn.load(file.path(R.home("lib"),"libRlapack.dylib"))
geigen <- function(A,B,jobvl=FALSE,jobvr=FALSE) {
# simplistic interface to Lapack dggev
# for generalized eigenvalue problem
# general matrices
# for symmetric matrices use dsygv
if(!is.matrix(A)) stop("Argument A should be a matrix")
if(!is.matrix(B)) stop("Argument B should be a matrix")
dimA <- dim(A)
if(dimA[1]!=dimA[2]) stop("A must be a square matrix")
dimB <- dim(B)
if(dimB[1]!=dimB[2]) stop("B must be a square matrix")
if(dimA[1]!=dimB[1]) stop("A and B must have the same dimensions")
if( is.complex(A) ) stop("A may not be complex")
if( is.complex(B) ) stop("B may not be complex")
jobvl.char <- if(jobvl) "V" else "N"
jobvr.char <- if(jobvr) "V" else "N"
n <- dimA[1]
# minimum amount of work memory
# for performance this needs to be set properly
# (see source dggev)
lwork <- as.integer(8*n)
work <- numeric(lwork)
if( jobvl && jobvr )
z <- .Fortran("dggev", jobvl.char, jobvr.char, n, A, n, B, n,
alphar=numeric(n), alphai=numeric(n),
beta=numeric(n), vl=matrix(0,nrow=n,ncol=n), n,
vr=matrix(0,nrow=n,ncol=n), n,
work, lwork,info=integer(1L))
else if(jobvl && !jobvr )
z <- .Fortran("dggev", jobvl.char, jobvr.char, n, A, n, B, n,
alphar=numeric(n), alphai=numeric(n),
beta=numeric(n), vl=matrix(0,nrow=n,ncol=n), n,
numeric(1), 1L,
work, lwork,info=integer(1L))
else if(!jobvl && jobvr )
z <- .Fortran("dggev", jobvl.char, jobvr.char, n, A, n, B, n,
alphar=numeric(n), alphai=numeric(n),
beta=numeric(n), numeric(1),1L,
vr=matrix(0,nrow=n,ncol=n), n,
work, lwork,info=integer(1L))
else
z <- .Fortran("dggev", jobvl.char, jobvr.char, n, A, n, B, n,
alphar=numeric(n), alphai=numeric(n),
beta=numeric(n), numeric(1), 1L, numeric(1), 1L,
work, lwork,info=integer(1L))
if( z$info > 0 ) stop(paste("Lapack dggev fails with info=",z$info))
# simplistic calculation of eigenvalues (see caveat in source dggev)
if( all(z$alphai==0) )
values <- z$alphar/z$beta
else
values <- complex(real=z$alphar, imaginary=z$alphai)/z$beta
if( jobvl && jobvr )
return(list(values=values, vl=z$vl, vr=z$vr))
else if( jobvl )
return(list(values=values, vl=z$vl))
else if( jobvr )
return(list(values=values, vr=z$vr))
else
return(list(values=values))
}
A <- matrix(c(1457.738, 1053.181, 1256.953,
1053.181, 1213.728, 1302.838,
1256.953, 1302.838, 1428.269), nrow=3, byrow=TRUE)
B <- matrix(c(4806.033, 1767.480, 2622.744,
1767.480, 3353.603, 3259.680,
2622.744, 3259.680, 3476.790), nrow=3, byrow=TRUE)
A
B
geigen(A, B)
geigen(A, B, TRUE , TRUE)
geigen(A, B, TRUE , FALSE)
geigen(A, B, FALSE, TRUE)
<end R code>
Berend
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
