## Wednesday, September 24, 2008

### [TECH] Computing Eigen Values and Eigen Vectors using LAPACK

I have used the LAPACK to compute the eigen values, till recently I needed the eigen vectors both left and right. I thought this could be really handy if you need to compute eigen values and eigen vectors. Download both library and this file from here

/*Computing both the eigen values and eigen vectors using LAPACK.
*vamsi.krishnak (at) gmail (dot) com
**/
#include "f2c.h"
#include "clapack.h"
#include<stdio.h>
#include<math.h>
#include<malloc.h>
#include<assert.h>

int main(){
integer dim=0;
unsigned int i;
doublereal *A = NULL;
real *eig_values = NULL;
real *work;
char JOBVL='V'; /*Left eigen vectors*/
char JOBVR='V'; /*No right eigen vectors*/
integer iwork[1];
integer lwork,liwork=1,info,lda;
doublereal *WR=NULL;
doublereal *WI=NULL;
doublereal *VL=NULL;
doublereal *VR=NULL;
integer LDVR=1;
integer LDVL=1;
integer LDA;
doublereal *WORK=NULL;
integer LWORK;

printf("Please Enter the dimension of the array:\n");
scanf("%d",&dim); LDA=dim; lda=dim; LDVL=dim;LDVR=dim;
eig_values = (real *) malloc(sizeof(real)*dim);
lwork = 10*dim;
work = (real *) malloc(sizeof(real)*lwork);
A = (doublereal *) malloc(sizeof(doublereal)*(dim*dim));
printf("Please Enter the elements of the matrix\n");
i = 0;

/*Double real and imaginary parts*/
WR = (doublereal *)malloc(sizeof(doublereal)*dim);
WI = (doublereal *)malloc(sizeof(doublereal)*dim);
VL = (doublereal *)malloc(sizeof(doublereal)*(dim*dim));
VR = (doublereal *)malloc(sizeof(doublereal)*(dim*dim));
LWORK = 6*dim;
WORK = (doublereal *)malloc(sizeof(doublereal)*LWORK);
assert(WI && WR && VL && VR);
/**************/
do{
scanf("%lf",&A[i]);
i++;
}while(i<(dim*dim));
/*Compute the Eigen values*/
//ssyev_("N","U",&dim,A,&dim,eig_values,work,&lwork/*,iwork,&liworki*/,&info);
dgeev_(&JOBVL,&JOBVR,&dim,A,&LDA,WR,WI,VL,&LDVL,VR,&LDVR,WORK,&LWORK,&info);
if(!info){
for(i=0;i<dim;i++){
printf("%lf+i%lf\n",WR[i],WI[i]);
}
printf("===Right Eigen Vectors===\n");
int k; int conj;
for(i=0;i<dim;i++){
/*Print i^th right eigen vector*/
if(fabs(WI[i] > 0)){
/*eigen vector i and eigen vector i+1*/
for(conj=0;conj<2;conj++){
printf("[ ");
for(k=0;k<dim;k++){
if(!conj){
printf("(%lf)+i(%lf) ",VL[i*dim+k],VL[(i+1)*dim+k]);

}else{
printf("(%lf)-i(%lf) ",VL[i*dim+k],VL[(i+1)*dim+k]);
}
}
printf("  ]\n");
}
i++;
}else{
/*real eigen vector*/
printf("[  ");
for(k=0;k<dim;k++){
printf("%lf ",VL[i*dim+k]);
}
printf("] \n");
}
printf("\n");
}
}
}