CODEWITHSUNDEEP
rlapack
rocastocks
rocastocks

Reputation: 153

Get eigenvalues from tridiagonal symmetric matrix in R like in lapack dstev

I was wondering if it is possible to diagonalize a symmetric tridiagonal matrix (real) in R passing only the diagonal vectors so you can avoid generating the complete matrix and use an optimized routine like dstev or dsteqr from Lapack

By the moment I know that for the R function eigen you can just specify whether the matrix is symmetric or not.

Is it worth to try to access directly to lapack or just do the a manual algorithm for this special case, in order to get good performance?

Upvotes: 2

Views: 576

Answers (1)

rocastocks
rocastocks

Reputation: 153

After some research and trial and error I found solutions to this problem by myself.

There are native ways to call external code in R, the routines .Fortran and .C can execute compiled code from those languages without the need of an intermediate layer.

The following code explains how to use the LAPACK routine dstev in order to calculate eigenvalues and eigenvectors from a symmetric matrix. To do this, it is necessary to have a BLAS/LAPACK implementation installed on your machine, such as Atlas, OpenBLAS or MKL.

dyn.load("/usr/lib/liblapack.so")  #Load Lapack
Jobz <- as.character('V')            #Get The eigenvalues and Eigenvectors
N <- as.integer(10)                  #Matrix Dimension
D <- runif(10)                       #Diagonal
E <- runif(9)                        #Upper (and lower) diagonal
Z <- array(0.0,N*N)                  #storage for eigenvectors
Ldz <- N                             #Leading dimension of Z (in our case number of rows)
Work <- array(0.0,2*N)               #Storage for processing
info <- as.integer(0)                #Info flag

#Call the fortran routine
res<-.Fortran("dstev", Jobz, N, D, E, Z, Ldz, Work, info)

#Extract the information
evalues <- res[[3]]
evectors <- matrix(res[[5]],N,N)

Upvotes: 2

Related Questions