# GSL sample: Eigenvalues and Engenvectors for Real Symmetric Matrix

Last Update: Oct. 4, 2007

### Sample code

eigen_symmetric.c

```/***********************************************/
/* 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;
}```

### Example execution

```\$ ./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 ```

<- Go back