#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
 
int main (void)
{
size_t 	i,j,NMax=1000;
size_t	EvalOut= 10;
double 	eval_i;
double 	*data=malloc(NMax*NMax*sizeof(double));
 
for(i=0; i<NMax; i++)
  {
  for(j=0; j<NMax; j++)
    {
    data[i*NMax + j] = 1.0/(i+j+1.0);
    }
  }   
 
 gsl_matrix_view m 
 = gsl_matrix_view_array (data, NMax, NMax);
 
 gsl_vector *eval = gsl_vector_alloc(NMax);
 gsl_matrix *evec = gsl_matrix_alloc(NMax, NMax);
 
 gsl_eigen_symmv_workspace * w = 
 gsl_eigen_symmv_alloc (NMax);
 
 gsl_eigen_symmv (&m.matrix, eval, evec, w);
 
 gsl_eigen_symmv_free (w);
 
 gsl_eigen_symmv_sort (eval, evec, 
 GSL_EIGEN_SORT_ABS_ASC);
 
 for(i=NMax-1; i>(NMax-EvalOut-1); i--)
 {
 eval_i = gsl_vector_get (eval, i);
 printf ("eigenvalue = %g\n", eval_i);
 }
 
 gsl_vector_free (eval);
 gsl_matrix_free (evec);
 
 return 0;
 }