gen_eigen.c

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

void 
general_eigen (double * A, double * B, double * x, double * lambda, int * size)
{
       size_t N = *size;
       gsl_matrix_view m 
         = gsl_matrix_view_array (A, N, N);

       gsl_matrix_view b 
         = gsl_matrix_view_array (B, N, N);
    
       gsl_vector_view eval
		 = gsl_vector_view_array (lambda, N);

       gsl_matrix_view evec
         = gsl_matrix_view_array (x, N, N);
    
       gsl_eigen_gensymmv_workspace * w = 
         gsl_eigen_gensymmv_alloc (N);
       
       gsl_eigen_gensymmv (&m.matrix, &b.matrix, &eval.vector, &evec.matrix, w);
     
       gsl_eigen_gensymmv_free (w);
     
       gsl_eigen_symmv_sort (&eval.vector, &evec.matrix, 
                             GSL_EIGEN_SORT_ABS_ASC);
            
}


int
main (void)
{
   double A[] = {		 1.0  , 1/2.0, 1/3.0, 1/4.0,
                         1/2.0, 1/3.0, 1/4.0, 1/5.0,
                         1/3.0, 1/4.0, 1/5.0, 1/6.0,
                         1/4.0, 1/5.0, 1/6.0, 1/7.0 };
   double B[] = {		 1.0  , 1/2.0, 1/3.0, 1/4.0,
                         1/2.0, 1/3.0, 1/4.0, 1/5.0,
                         1/3.0, 1/4.0, 1/5.0, 1/6.0,
                         1/4.0, 1/5.0, 1/6.0, 1/7.0 };
	double x[16];
	double lambda[4];
	int n = 4;
	general_eigen(A,B, x, lambda, &n);

	return 0;
}

general_eigen.R

# example

dyn.load("gen_eigen.so")
general_eigen <- function(A, B){
	n <- sqrt(length(A))
	out <- .C("general_eigen", as.double(A), as.double(B), double(n^2), double(n), as.integer(n))
	list(evals = out[[4]], evecs = out[[3]])
}


  A <- c( 0.24,   0.39 ,   0.42 ,  -0.16 ,  
           0.39 ,  -0.11 ,   0.79 ,   0.63 ,  
           0.42 ,   0.79 ,  -0.25 ,   0.48 ,  
          -0.16 ,   0.63 ,   0.48 ,  -0.03 )

B <- c(     4.16 ,  -3.12 ,   0.56 ,  -0.10 ,  
          -3.12 ,   5.03 ,  -0.83 ,   1.09 ,  
           0.56 ,  -0.83 ,   0.76 ,   0.34 ,  
          -0.10 ,   1.09 ,   0.34 ,   1.18)
o

Makevars

PKG_CPPFLAGS=-Wall 
PKG_LIBS=-lm -lgsl -lgslcblas 

Generalized Eigenvalue problem in R

A simple R wrapper to the gsl generalized eigenvalue problem, http://www.gnu.org/software/gsl/manual/html_node/Real-Generalized-Symmetric_002dDefinite-Eigensystems.html

Requires the gsl library. Makevars file specifies the linkage to the gsl libraries, then the R CMD SHLIB command will generate the shared object library called by the R code.