/*--------------------------------------------------------------------
  This program computes the eign values and eign vectors of matrix
  using the GSL libraries.
  
  The matrix for which the eign values & eign vectors are calculated,
  has to be given at the run time in a '.mat' file.

  This program has a test part also which makes sure that you get the 
  correct answers.
  
  
  Jayanti Prasad, April 27, 2011
  
  Commments & Suggestion : prasad.jayanti@gmail.com		
  
  -------------------------------------------------------------------- */

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

int test_eign(int , double [], double [], double[]); 

int  main (int argc, char *argv[])    {
  double *data,*data1; 
  FILE *inpf;
  int i, j,ndim; 
 
  if (argc < 2){
    fprintf(stderr,"./compute_inverse <matrix file>\n");
    return(-1);  
  }
 
  inpf=fopen(argv[1],"r");
  fscanf(inpf,"%d\n",&ndim);
  data=(double *)malloc(ndim*ndim*sizeof(double));
  data1=(double *)malloc(ndim*ndim*sizeof(double));

  for(i=0; i < ndim; i++){
    for(j=0; j < ndim; j++){
      fscanf(inpf,"%lf",&data[j+ndim*i]);
        data1[j+ndim*i]=data[j+ndim*i];
       }
    fscanf(inpf,"\n");
    
  }
  fclose(inpf);
 
 
  gsl_matrix_view m = gsl_matrix_view_array (data, ndim, ndim);
  // fill the matrix for which eign values/vectors are computed
  
  gsl_vector *eval  = gsl_vector_alloc (ndim);
  // a vector which contains the eign values 

  gsl_matrix *evec  = gsl_matrix_alloc (ndim, ndim);
  // a mtrix colums of which are eign vectors 
  
  gsl_eigen_symmv_workspace * w =  gsl_eigen_symmv_alloc (ndim);
  // this is just for work space 
  
  gsl_eigen_symmv (&m.matrix, eval, evec, w);
  // pass the matrix to main engine and get back the 
  // eign values and eign vectors 
  
  gsl_eigen_symmv_free (w);
  
  gsl_eigen_symmv_sort (eval, evec,GSL_EIGEN_SORT_ABS_ASC);
  
  // this section writies out put i.e., eign values & vectors, one by one 
  

  printf(" ----- Eign values \n");
  for(i=0; i < ndim ; i++)
    printf(" %2.6f ",gsl_vector_get(eval,i));
  printf("\n ---- Eign Vectors (Columns) \n");

  for(i=0; i < ndim; i++){
    for(j=0; j < ndim; j++)
      printf(" %2.6f ",gsl_matrix_get(evec,i,j));
    printf("\n");
  
  }
  printf("\n");
 

  { 
  double P[ndim*ndim],lambda[ndim];

  for(i=0; i < ndim; i++)
     lambda[i] = gsl_vector_get(eval,i); 
  for(i=0; i < ndim; i++)
    for(j=0; j < ndim; j++)
    P[j+ndim*i] = gsl_matrix_get(evec,i,j); 
 
  test_eign(ndim,lambda,data1,P);
}

  gsl_vector_free (eval);
  gsl_matrix_free (evec);
 
  return 0;
}

int test_eign(int ndim, double l[], double A[], double  P[]){
  int i, j, k; 
  double  *I,*B,*X,*Y;
  I=(double *)malloc(ndim*ndim*sizeof(double));
  B=(double *)malloc(ndim*ndim*sizeof(double));
  Y=(double *)malloc(ndim*sizeof(double));
  X=(double *)malloc(ndim*sizeof(double));

 printf("-------Checking -----------------------------\n"); 

  for(i=0; i < ndim; i++){
   for(j=0; j < ndim; j++){
     if (i==j)
        I[j+ndim*i]=1.0;
     else
        I[j+ndim*i]=0.0; 
   }
  }

  for(i=0; i < ndim; i++){ 
      for(j=0; j < ndim ; j++){
        for(k=0; k < ndim ; k++){
          B[k+ndim*j]=A[k+ndim*j]-I[k+ndim*j] *l[i];
       }
    }

   for(j=0; j < ndim ; j++)
       X[j]=P[i+ndim*j];


   for(j=0; j < ndim; j++){
       Y[j] = 0.0;
      for(k=0; k < ndim; k++){
        Y[j] += B[k+ndim*j] * X[k];       
    }
   }

  printf("Eign Value=%2.6f\n",l[i]); 
  printf("Eign Vector =  (" );
  for(j=0; j < ndim ; j++)
       printf("  %2.6f ",X[j]);
   printf(")\n" );

  printf("Residual:  (A-lambda *I)*X = (" );
  for(j=0; j < ndim ; j++)
       printf("  %2.6f ",Y[j]);
   printf(")\n" );
}

  return(0); 
}