#include <stdio.h>
     #include <gsl/gsl_math.h>
     #include <gsl/gsl_eigen.h>
     
     int
     main (void)
     {
       double data[] = { -1.0, 1.0, -1.0, 1.0,
                         -8.0, 4.0, -2.0, 1.0,
                         27.0, 9.0, 3.0, 1.0,
                         64.0, 16.0, 4.0, 1.0 };
     
       gsl_matrix_view m 
         = gsl_matrix_view_array (data, 4, 4);
     
       gsl_vector_complex *eval = gsl_vector_complex_alloc (4);
       gsl_matrix_complex *evec = gsl_matrix_complex_alloc (4, 4);
     
       gsl_eigen_nonsymmv_workspace * w = 
         gsl_eigen_nonsymmv_alloc (4);
       
       gsl_eigen_nonsymmv (&m.matrix, eval, evec, w);
     
       gsl_eigen_nonsymmv_free (w);
     
       gsl_eigen_nonsymmv_sort (eval, evec, 
                                GSL_EIGEN_SORT_ABS_DESC);
       
       {
         int i, j;
     
         for (i = 0; i < 4; i++)
           {
             gsl_complex eval_i 
                = gsl_vector_complex_get (eval, i);
             gsl_vector_complex_view evec_i 
                = gsl_matrix_complex_column (evec, i);
     
             printf ("eigenvalue = %g + %gi\n",
                     GSL_REAL(eval_i), GSL_IMAG(eval_i));
             printf ("eigenvector = \n");
             for (j = 0; j < 4; ++j)
               {
                 gsl_complex z = 
                   gsl_vector_complex_get(&evec_i.vector, j);
                 printf("%g + %gi\n", GSL_REAL(z), GSL_IMAG(z));
               }
           }
       }
     
       gsl_vector_complex_free(eval);
       gsl_matrix_complex_free(evec);
     
       return 0;
     }