| #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; |
| } |