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

[Berend Hasselman]

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>
> 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