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 <imad@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 <bhh@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: jgrn307@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
http://www.geog.illinois.edu/people/JonathanGreenberg.html

	[[alternative HTML version deleted]]

