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