htswanalysis/MACS/lib/gsl/gsl-1.11/doc/examples/eigen.c

   1 #include <stdio.h>
   2 #include <gsl/gsl_math.h>
   3 #include <gsl/gsl_eigen.h>
   4
   5 int
   6 main (void)
   7 {
   8   double data[] = { 1.0  , 1/2.0, 1/3.0, 1/4.0,
   9                     1/2.0, 1/3.0, 1/4.0, 1/5.0,
  10                     1/3.0, 1/4.0, 1/5.0, 1/6.0,
  11                     1/4.0, 1/5.0, 1/6.0, 1/7.0 };
  12
  13   gsl_matrix_view m
  14     = gsl_matrix_view_array (data, 4, 4);
  15
  16   gsl_vector *eval = gsl_vector_alloc (4);
  17   gsl_matrix *evec = gsl_matrix_alloc (4, 4);
  18
  19   gsl_eigen_symmv_workspace * w =
  20     gsl_eigen_symmv_alloc (4);
  21
  22   gsl_eigen_symmv (&m.matrix, eval, evec, w);
  23
  24   gsl_eigen_symmv_free (w);
  25
  26   gsl_eigen_symmv_sort (eval, evec,
  27                         GSL_EIGEN_SORT_ABS_ASC);
  28
  29   {
  30     int i;
  31
  32     for (i = 0; i < 4; i++)
  33       {
  34         double eval_i
  35            = gsl_vector_get (eval, i);
  36         gsl_vector_view evec_i
  37            = gsl_matrix_column (evec, i);
  38
  39         printf ("eigenvalue = %g\n", eval_i);
  40         printf ("eigenvector = \n");
  41         gsl_vector_fprintf (stdout,
  42                             &evec_i.vector, "%g");
  43       }
  44   }
  45
  46   gsl_vector_free (eval);
  47   gsl_matrix_free (evec);
  48
  49   return 0;
  50 }