1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
|
#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;
}
|
