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

}

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

Subscribe to:
Post Comments (Atom)

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.

ReplyDeleteThanks

The Makefile is something like:

ReplyDelete#

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

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

ReplyDeleteMy example uses zge functions.

z: double complex

ge: general matrices.

hi,

ReplyDeletei 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

Hi,

ReplyDeleteYour 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

Hi!! Thanks!! nice example!!

ReplyDeleteThank you very much. This is quite useful.

ReplyDelete