Last Update: Oct. 4, 2007
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 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 | /***********************************************/ /* T.Kouya's GSL sample program collection */ /* Eingenvalue Problem */ /* for Real Symmetrix Matrix */ /* */ /* Written by Tomonori Kouya */ /* */ /* Version 0.01: 2007-10-04 */ /***********************************************/ #include <stdio.h> #include <gsl/gsl_math.h> #include <gsl/gsl_blas.h> #include <gsl/gsl_vector.h> #include <gsl/gsl_matrix.h> #include <gsl/gsl_eigen.h> #include <gsl/gsl_sort_vector.h> /* Dimension of Matrix and Vectors */ #define DIM 5 int main( void ) { int i, j; gsl_permutation *perm; gsl_matrix *frank_a, *matrix, *eig_vecs, *tmp_matrix; gsl_vector *eig_vec, *eig; gsl_eigen_symm_workspace *workspace; gsl_eigen_symmv_workspace *workspace_v; /* allocate a, x, b */ frank_a = gsl_matrix_alloc (DIM, DIM); matrix = gsl_matrix_alloc (DIM, DIM); eig_vec = gsl_vector_alloc (DIM); eig = gsl_vector_alloc (DIM); /* set Frank matrix */ /* a_ij = DIM - min(i,j) + 1 */ for (i = 0; i < DIM; i++) for (j = 0; j < DIM; j++) gsl_matrix_set (frank_a, i, j, ( double )DIM - ( double ) GSL_MAX (i, j) ); /* Print matrix */ printf ( "Frank Matrix(DIM = %d)\n" , DIM); for (i = 0; i < DIM; i++) { printf ( "%3d: " , i); for (j = 0; j < DIM; j++) printf ( "%g " , gsl_matrix_get (frank_a, i, j)); printf ( "\n" ); } printf ( "\n" ); /* Get Eigensystem of Frank Matrix : Eigenvalues ONLY */ /* Initialize workspace and copy Frank matrix to "matrix" */ workspace = gsl_eigen_symm_alloc (DIM); gsl_matrix_memcpy (matrix, frank_a); /* get eigenvalues (matrix is destroyed) */ gsl_eigen_symm (matrix, eig, workspace); /* sort eigenvalues */ gsl_sort_vector (eig); /* print eigenvalues */ printf ( " i eigenvalues \n" ); for (i = 0; i < DIM; i++) printf ( "%3d %25.17e\n" , i, gsl_vector_get (eig, i)); /* free workspace */ gsl_eigen_symm_free (workspace); /* Get Eigensystem of Frank Matrix : Eigenvalues and Eigenvectors */ /* Initialize workspace and eigenvectors stored, and set matrix */ eig_vecs = gsl_matrix_alloc (DIM, DIM); workspace_v = gsl_eigen_symmv_alloc (DIM); gsl_matrix_memcpy (matrix, frank_a); /* get eigenvalues and its eigenvectors (matrix is destroyed) */ gsl_eigen_symmv (matrix, eig, eig_vecs, workspace_v); /* sort eigenvalues with permutation */ perm = gsl_permutation_alloc (DIM); gsl_sort_vector_index (perm, eig); /* print eigenvalues */ printf ( " i eigenvalues \n" ); for (i = 0; i < DIM; i++) printf ( "%3d %25.17e\n" , i, gsl_vector_get (eig, gsl_permutation_get (perm, i))); printf ( "\n" ); /* print eigenvectors */ printf ( "eigenvectors : [v1 v2 ... vn]\n" ); for (i = 0; i < DIM; i++) printf ( "eig%2d = %5.1e " , i, gsl_vector_get (eig, i)); printf ( "\n" ); for (i = 0; i < DIM; i++) { for (j = 0; j < DIM; j++) printf ( "%16.9e " , gsl_matrix_get (eig_vecs, i, j)); printf ( "\n" ); } printf ( "\n" ); /* free permutation */ gsl_permutation_free (perm); /* A * [v1 v2 ... vn] - diag[eig1 eig2 ... eign] * [v1 v2 ... vn] == 0 ? */ tmp_matrix = gsl_matrix_alloc (DIM, DIM); gsl_matrix_set_identity (tmp_matrix); for (i = 0; i < DIM; i++) for (j = 0; j < DIM; j++) gsl_matrix_set (matrix, j, i, gsl_vector_get (eig, i)); gsl_matrix_mul_elements (matrix, eig_vecs); gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, frank_a, eig_vecs, -1.0, matrix); /* print residual */ printf ( "A * [v1 v2 ... vn] - diag[eig1 eig2 ... eign] * [v1 v2 ... vn]\n" ); for (i = 0; i < DIM; i++) { printf ( "%3d: " , i); for (j = 0; j < DIM; j++) printf ( "%g " , gsl_matrix_get (matrix, i, j)); printf ( "\n" ); } printf ( "\n" ); /* free workspace */ gsl_matrix_free (eig_vecs); gsl_matrix_free (tmp_matrix); gsl_eigen_symmv_free (workspace_v); /* free a, x, b */ gsl_matrix_free (frank_a); gsl_matrix_free (matrix); gsl_vector_free (eig_vec); gsl_vector_free (eig); return 0; } |
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 | $ ./eigen_symmetric Frank Matrix(DIM = 5) 0: 5 4 3 2 1 1: 4 4 3 2 1 2: 3 3 3 2 1 3: 2 2 2 2 1 4: 1 1 1 1 1 i eigenvalues 0 2.71554129338822337e-01 1 3.53253282893739418e-01 2 5.82964498293740307e-01 3 1.44869056979664323e+00 4 1.23435375196770654e+01 i eigenvalues 0 2.71554129338822337e-01 1 3.53253282893739418e-01 2 5.82964498293740307e-01 3 1.44869056979664323e+00 4 1.23435375196770654e+01 eigenvectors : [v1 v2 ... vn] eig 0 = 1.2e+01 eig 1 = 1.4e+00 eig 2 = 5.8e-01 eig 3 = 3.5e-01 eig 4 = 2.7e-01 5.968847877e-01 5.485287320e-01 -4.557341407e-01 3.260186796e-01 1.698911240e-01 5.485287320e-01 1.698911240e-01 3.260186796e-01 -5.968847877e-01 -4.557341407e-01 4.557341407e-01 -3.260186796e-01 5.485287320e-01 1.698911240e-01 5.968847877e-01 3.260186796e-01 -5.968847877e-01 -1.698911240e-01 4.557341407e-01 -5.485287320e-01 1.698911240e-01 -4.557341407e-01 -5.968847877e-01 -5.485287320e-01 3.260186796e-01 A * [v1 v2 ... vn] - diag[eig1 eig2 ... eign] * [v1 v2 ... vn] 0: -3.13638e-15 -3.33067e-16 -7.77156e-16 -6.66134e-16 4.996e-16 1: -5.80092e-15 -3.33067e-16 -1.66533e-15 1.11022e-15 2.27596e-15 2: -4.02456e-15 5.55112e-16 -9.99201e-16 -4.44089e-16 5.55112e-17 3: -1.02696e-15 3.33067e-16 -8.88178e-16 2.22045e-16 8.32667e-16 4: -1.36002e-15 -1.11022e-16 -8.88178e-16 -1.11022e-16 3.88578e-16 |