Sunday, May 2, 2010

Eigenvalues: C/Lapack

Here is another example on the use of Lapack. This time the objective is to calculate the eigenvectors and eigenvalues of a complex matrix

#include<stdio.h>
#include<math.h>

#include<complex.h>
#include <stdlib.h>

//.......................................................................................................
void zgeTranspose( double complex *Transposed, double complex *M ,int n)
{


int
i,j;
for
(i=0;i<n;i++)

for
(j=0;j<n;j++) Transposed[i+n*j] = M[i*n+j];
}


//......................................................................................................
// MatrixComplexEigensystem: computes the eigenvectors and eigenValues of input matrix A
// The eigenvectors are stored in columns
//.....................................................................................................
void MatrixComplexEigensystem( double complex *eigenvectorsVR, double complex *eigenvaluesW, double complex *A, int N)
{


int
i;

double
complex *AT = (double complex*) malloc( N*N*sizeof(double complex) );

zgeTranspose( AT, A , N);

char
JOBVL ='N'; // Compute Right eigenvectors

char JOBVR ='V'; // Do not compute Left eigenvectors

double
complex VL[1];

int
LDVL = 1;
int
LDVR = N;

int
LWORK = 4*N;

double
complex *WORK = (double complex*)malloc( LWORK*sizeof(double complex));

double
complex *RWORK = (double complex*)malloc( 2*N*sizeof(double complex));


int
INFO;

zgeev_( &JOBVL, &JOBVR, &N, AT , &N , eigenvaluesW ,

VL, &LDVL,
eigenvectorsVR, &LDVR,
WORK,
&
LWORK, RWORK, &INFO );

zgeTranspose( AT, eigenvectorsVR , N);

for
(i=0;i<N*N;i++) eigenvectorsVR[i]=AT[i];


free(WORK);
free(RWORK);
free(AT);
}



int
main()
{

int
i,j;
const
int N = 3;


double
complex A[] = { 1.+I , 2. , 3 , 4. , 5.+I , 6. , 7., 8., 9. + I};

double
complex eigenVectors[N*N];
double
complex eigenValues[N];


MatrixComplexEigensystem( eigenVectors, eigenValues, A, N);

printf("\nEigenvectors\n");

for
(i=0;i<N;i++){
for
(j=0;j<N;j++) printf(" (%f,%f) \t", eigenVectors[i*N + j]);

printf("\n");
}


printf("\nEigenvalues \n");
for
(i=0;i<N;i++) printf("\n (%f, %f) \t", eigenValues[i] );


printf("\n------------------------------------------------------------\n");
return
0;

}

7 comments:

  1. A very useful example code that I have been searching for days. Could you also please add a few lines on how to compile this code. Also, are there any specialised functions to diagonalize real symmetric matrices using Lapack in C.
    Thanks

    ReplyDelete
  2. The Makefile is something like:

    #
    CC = gcc

    #edit LAPACK_PATH if necessary
    LAPACK_PATH = /usr/lib64/atlas

    a,out: main.c
    $(CC) main.c -L$(LAPACK_PATH) -llapack -lblas -lgfortran -lm

    ReplyDelete
  3. Yes, Lapack contains methods for special kinds of matrices in order to optimize the code. You can find the function you want for example at http://www.netlib.org/lapack/double/.

    My example uses zge functions.
    z: double complex
    ge: general matrices.

    ReplyDelete
  4. hi,
    i am trying to compile this code but when i am doing
    gcc main.c -llapack -lblas -lgfortran -lm
    i am getting an error as
    /usr/bin/ld: cannot find -llapack
    Please if you can guide me in removing this error

    Thanks
    Rahi

    ReplyDelete
  5. Hi,

    Your compiler cannot find the lapack library file
    which can be
    liblapack.so or liblapack.a

    make sure you have those files.

    You may need to edit
    /etc/ld.so.conf
    to include the path where the libraries are located
    and execute ldconfig

    ReplyDelete
  6. Thank you very much. This is quite useful.

    ReplyDelete