/*    ATK IV - Numerical Programming (2006)      */
/*    Excercise 9.1 -  Eigenvalue Problems       */

#include<stdio.h>
#include<math.h>
#define NRANSI
#include<nr.h>
#include<nrutil.h>

int main(void){
  float **a, **v;
  float *d;
  int n,nrot;
  float iprod=0.0;

  FILE *file;
  int i,j;

  file=fopen("matrix_sym.dat","rt");
  fscanf(file,"%d\n",&n);

  a=matrix(1,n,1,n);
  v=matrix(1,n,1,n);
  d=vector(1,n);

  for(i=1;i<=n;i++){
    for(j=1;j<=n;j++){
      fscanf(file,"%f",&a[i][j]);
    }    
  }

  jacobi(a,n,d,v,&nrot);
  eigsrt(d,v,n);
  
  printf("\n");

  for(i=1;i<=n;i++){
    printf("Eigenvalue: %f\n",d[i]);
    printf("Corresponding eigenvector:  ");
    for(j=1;j<=n;j++){
      printf("%f ",v[j][i]);
    }
    printf("\n");
  }

  for(i=1;i<=n;i++){
    iprod += v[i][1]*v[i][2];
  }
  printf("\nInner product between 1st and 2nd eigenvector: %f\n",iprod);
  free_matrix(a,1,n,1,n);
  free_matrix(v,1,n,1,n);
  free_vector(d,1,n);
  return 0;

}