Building in Berend's suggestions I think this function should work for most people (I'm going to wrap it into a package but figured people may want to grab this directly):
# Please see http://www.netlib.org/lapack/double/dggev.f for a description of inputs and outputs. Rdggev <- function(A,B,JOBVL=FALSE,JOBVR=TRUE) { # 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 <i...@estarcion.net> # Contributions from Berend Hasselman. if( .Platform$OS.type == "windows" ) { Lapack.so <- file.path(R.home("bin"),paste("Rlapack",.Platform$dynlib.ext,sep="")) } else { Lapack.so <- file.path(R.home("modules"),paste("lapack",.Platform$dynlib.ext,sep="")) } dyn.load(Lapack.so) if(JOBVL) { JOBVL="V" } else { JOBVL="N" } if(JOBVR) { JOBVR="V" } else { JOBVR="N" } 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") # 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") # simplistic calculation of eigenvalues (see caveat in http://www.netlib.org/lapack/double/dggev.f) if( all(Rdggev_out$ALPHAI==0) ) Rdggev_out$GENEIGENVALUES <- Rdggev_out$ALPHAR/Rdggev_out$BETA else Rdggev_out$GENEIGENVALUES <- complex(real=Rdggev_out$ALPHAR, imaginary=Rdggev_out$ALPHAI)/Rdggev_out$BETA return(Rdggev_out) } Thanks! --j On Sun, Apr 22, 2012 at 2:27 PM, Berend Hasselman <b...@xs4all.nl> wrote: > > On 22-04-2012, at 21:08, Jonathan Greenberg wrote: > > > Thanks all (particularly to you, Berend) -- I'll push forward with these > solutions and integrate them into my code. I did come across geigen while > rooting around in the CCA code but its not formally documented (it just > says "for internal use" or something along those lines) and as you found > out above, it does not produce the same solution as the dggev. It would be > nice to have a more complete set of formal packages for doing LA in R > (rather than having to hand-write .Fortran calls) but I'll leave that to > someone with more expertise in linear algebra than me. Something that > perhaps matches the SciPy set of functions (both in terms of input and > output): > > > > http://docs.scipy.org/doc/scipy/reference/linalg.html > > > > Some of these are already implemented, but clearly not all of them. > > Package CCA has package fda as dependency. > And package fda defines a function geigen. > The first 14 lines of this function are > > geigen <- function(Amat, Bmat, Cmat) > { > # solve the generalized eigenanalysis problem > # > # max {tr L'AM / sqrt[tr L'BL tr M'CM] w.r.t. L and M > # > # Arguments: > # AMAT ... p by q matrix > # BMAT ... order p symmetric positive definite matrix > # CMAT ... order q symmetric positive definite matrix > # Returns: > # VALUES ... vector of length s = min(p,q) of eigenvalues > # LMAT ... p by s matrix L > # MMAT ... q by s matrix M > > It's not clear to me how it is used and exactly what it is doing and how > that compares with Lapack. > > Berend > > -- Jonathan A. Greenberg, PhD Assistant Professor Department of Geography and Geographic Information Science University of Illinois at Urbana-Champaign 607 South Mathews Avenue, MC 150 Urbana, IL 61801 Phone: 415-763-5476 AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007 http://www.geog.illinois.edu/people/JonathanGreenberg.html [[alternative HTML version deleted]] ______________________________________________ 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.
