[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;
    printf("Please enter the Matrix A\n");
    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");
        }
    }
}