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

Paul Johnson pauljohn32 at gmail.com
Sat Apr 21 00:06:01 CEST 2012 

On Fri, Apr 20, 2012 at 2:18 PM, Jonathan Greenberg <jgrn at illinois.edu> 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


There's documentation on this in the "Writing R Extensions" manual,
See 1.2.1 Using Makevars and Section 6.7 "Numerical analysis
subrouties." They recommend looking at the package fastICA. It appears
to me you will use some platform neutral Makefile variables.

If you give a full working example from windows--your function, data
to use it--I'll see what I can do on the Linux side.  I've never had
need to worry about this question before, but it is about time I
learned.

Before you think about packaging something like this, I'd say step one
is to clean up your code so it is easier to read.  R CMD check won't
let you get bye with T or F in place of TRUE and FALSE, and you seem
to mix <- and = in your code, which makes it difficult to read.

pj


> "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)
> }
>
>
>
> On Fri, Apr 20, 2012 at 12:01 PM, Jonathan Greenberg <jgrn at illinois.edu>wrote:
>
>> Incidentally, if you want to test this out:
>>
>> > A
>>          [,1]     [,2]     [,3]
>> [1,] 1457.738 1053.181 1256.953
>> [2,] 1053.181 1213.728 1302.838
>> [3,] 1256.953 1302.838 1428.269
>>
>> > B
>>          [,1]     [,2]     [,3]
>> [1,] 4806.033 1767.480 2622.744
>> [2,] 1767.480 3353.603 3259.680
>> [3,] 2622.744 3259.680 3476.790
>>
>> I THINK this is correct for the other parameters:
>>         JOBVL="N"
>>         JOBVR="N"
>> N=dim(A)[[1]]
>>  LDA=N
>> LDB=N
>> LDVL=N
>>  LDVR=N
>> LWORK=max(1,8*N)
>>
>> .Call("La_dgeev", JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR=NA, ALPHAI=NA,
>> BETA=NA,
>>  VL=NA, LDVL, VR=NA, LDVR, WORK=NA, LWORK, INFO=NA, PACKAGE="base")
>>  LAPACK's documentation is:
>> http://www.netlib.org/lapack/double/dggev.f
>>
>> --j
>>
>> On Fri, Apr 20, 2012 at 11:58 AM, Jonathan Greenberg <jgrn at illinois.edu>wrote:
>>
>>> So I am a complete noob at doing this type of linking.  When I write this
>>> statement (all the input values are assigned, but I have to confess I don't
>>> know what to do with the output variables -- in this call they are all
>>> assigned "NA"):
>>>
>>> .Call("La_dgeev", JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR=NA, ALPHAI=NA,
>>> BETA=NA,
>>> + VL=NA, LDVL, VR=NA, LDVR, WORK=NA, LWORK, INFO=NA, PACKAGE="base")
>>>
>>> I get:
>>>
>>> Error in .Call("La_dgeev", JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR = NA,
>>>  :
>>>   C symbol name "La_dgeev" not in DLL for package "base"
>>>
>>> I'm running this on Windows 7 x64 version of R 2.14.2.  The
>>> R_ext/Lapack.h file states:
>>>
>>> /* DGEEV - compute for an N-by-N real nonsymmetric matrix A, the */
>>> /* eigenvalues and, optionally, the left and/or right eigenvectors */
>>> La_extern void
>>> F77_NAME(dgeev)(const char* jobvl, const char* jobvr,
>>> const int* n, double* a, const int* lda,
>>> double* wr, double* wi, double* vl, const int* ldvl,
>>>  double* vr, const int* ldvr,
>>> double* work, const int* lwork, int* info);
>>>
>>> Any ideas?
>>>
>>> --j
>>>
>>>
>>>
>>> On Fri, Apr 20, 2012 at 1:37 AM, Berend Hasselman <bhh at xs4all.nl> wrote:
>>>
>>>>
>>>> On 19-04-2012, at 20:50, Jonathan Greenberg wrote:
>>>>
>>>> > Folks:
>>>> >
>>>> > I'm trying to port some code from python over to R, and I'm running
>>>> into a
>>>> > wall finding R code that can solve a generalized eigenvalue problem
>>>> > following this function model:
>>>> >
>>>> >
>>>> http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eig.html
>>>> >
>>>> > Any ideas?  I don't want to call python from within R for various
>>>> reasons,
>>>> > I'd prefer a "native" R solution if one exists.  Thanks!
>>>>
>>>> An old thread is here:
>>>> http://tolstoy.newcastle.edu.au/R/help/05/06/6802.html
>>>>
>>>> R does provide eigen().
>>>> R has the Lapack routines dggev and dgges in its library.
>>>> You'd have to write your own .Fortran interface to those routines.
>>>>
>>>> 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 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
>>> http://www.geog.illinois.edu/people/JonathanGreenberg.html
>>>
>>
>>
>>
>> --
>> 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 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
>> http://www.geog.illinois.edu/people/JonathanGreenberg.html
>>
>
>
>
> --
> 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 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
> http://www.geog.illinois.edu/people/JonathanGreenberg.html
>
>        [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at 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.



-- 
Paul E. Johnson
Professor, Political Science    Assoc. Director
1541 Lilac Lane, Room 504     Center for Research Methods
University of Kansas               University of Kansas
http://pj.freefaculty.org            http://quant.ku.edu

More information about the R-help mailing list