htswanalysis/MACS/lib/gsl/gsl-1.11/doc/examples/eigen_nonsymm.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.0, -1.0, 1.0,
   9                     -8.0, 4.0, -2.0, 1.0,
  10                     27.0, 9.0, 3.0, 1.0,
  11                     64.0, 16.0, 4.0, 1.0 };
  12
  13   gsl_matrix_view m
  14     = gsl_matrix_view_array (data, 4, 4);
  15
  16   gsl_vector_complex *eval = gsl_vector_complex_alloc (4);
  17   gsl_matrix_complex *evec = gsl_matrix_complex_alloc (4, 4);
  18
  19   gsl_eigen_nonsymmv_workspace * w =
  20     gsl_eigen_nonsymmv_alloc (4);
  21
  22   gsl_eigen_nonsymmv (&m.matrix, eval, evec, w);
  23
  24   gsl_eigen_nonsymmv_free (w);
  25
  26   gsl_eigen_nonsymmv_sort (eval, evec,
  27                            GSL_EIGEN_SORT_ABS_DESC);
  28
  29   {
  30     int i, j;
  31
  32     for (i = 0; i < 4; i++)
  33       {
  34         gsl_complex eval_i
  35            = gsl_vector_complex_get (eval, i);
  36         gsl_vector_complex_view evec_i
  37            = gsl_matrix_complex_column (evec, i);
  38
  39         printf ("eigenvalue = %g + %gi\n",
  40                 GSL_REAL(eval_i), GSL_IMAG(eval_i));
  41         printf ("eigenvector = \n");
  42         for (j = 0; j < 4; ++j)
  43           {
  44             gsl_complex z = gsl_vector_complex_get(&evec_i.vector, j);
  45             printf("%g + %gi\n", GSL_REAL(z), GSL_IMAG(z));
  46           }
  47       }
  48   }
  49
  50   gsl_vector_complex_free(eval);
  51   gsl_matrix_complex_free(evec);
  52
  53   return 0;
  54 }