Last Update: Oct. 4, 2007
/***********************************************/
/* 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;
}
$ ./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