[R] Solve an ordinary or generalized eigenvalue problem in R?

Berend Hasselman bhh at xs4all.nl
Sat Apr 21 11:40:04 CEST 2012

```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 <jgrn at illinois.edu>
>
> 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.

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

```