# GSL sample: Eigenvalues and Engenvectors for Real Non-symmetric Matrix

Last Update: Aug. 29, 2008

### Sample code

eigen_nonsymmetric.c

```/***********************************************/
/* T.Kouya's GSL sample program collection     */
/* Eingenvalue Problem                         */
/*               for Real Non-symmetrix Matrix */
/*                                             */
/*                   Written by Tomonori Kouya */
/*                                             */
/* Version 0.01: 2008-08-29                    */
/***********************************************/
#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>

/* Dimension of Matrix and Vectors */
#define DIM 5

int main(void)
{
int i, j;
gsl_matrix *companion_a, *matrix;
gsl_matrix_complex *eig_vecs, *c_matrix, *tmp_matrix;
gsl_vector_complex *eig_vec, *eig;
gsl_eigen_nonsymm_workspace *workspace;
gsl_eigen_nonsymmv_workspace *workspace_v;

/* to extract elements from complex_vector & matrix */
gsl_complex c_element, alpha, beta;

/* allocate a, x, b */
companion_a = gsl_matrix_alloc(DIM, DIM);
matrix = gsl_matrix_alloc(DIM, DIM);
eig_vec = gsl_vector_complex_alloc(DIM);
eig = gsl_vector_complex_alloc(DIM);

/* set companion matrix: (DIM + 1) * x^DIM + ... + 2 * x + 1 = 0 */
/* a_1j = (DIM - j) / DIM */
/* a_i,i-1 = 1 */
for(i = 0; i < DIM; i++)
gsl_matrix_set(companion_a, 0, i, -(double)(DIM - i) / (DIM + 1));
for(i = 1; i < DIM; i++)
{
for(j = 0; j < DIM; j++)
gsl_matrix_set(companion_a, i, j, 0.0);
gsl_matrix_set(companion_a, i, i - 1, 1.0);
}

/* Print matrix */
printf("Companion 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(companion_a, i, j));
printf("\n");
}
printf("\n");

/* Get Eigensystem of Companion Matrix : Eigenvalues ONLY */

/* Initialize workspace and copy companion matrix to "matrix" */
workspace = gsl_eigen_nonsymm_alloc(DIM);
gsl_matrix_memcpy(matrix, companion_a);

/* get eigenvalues (matrix is destroyed) */
gsl_eigen_nonsymm(matrix, eig, workspace);

/* print eigenvalues */
printf(" i       eigenvalues        \n");
for(i = 0; i < DIM; i++)
{
c_element = gsl_vector_complex_get(eig, i);
printf("%3d %25.17e + %25.17e * i\n", i, GSL_REAL(c_element), GSL_IMAG(c_element));
}
printf("\n");

/* free workspace */
gsl_eigen_nonsymm_free(workspace);

/* Get Eigensystem of Companion Matrix : Eigenvalues and Eigenvectors */

/* Initialize workspace and eigenvectors stored, and set matrix */
eig_vecs = gsl_matrix_complex_alloc(DIM, DIM);
workspace_v = gsl_eigen_nonsymmv_alloc(DIM);
gsl_matrix_memcpy(matrix, companion_a);

/* get eigenvalues and its eigenvectors (matrix is destroyed) */
gsl_eigen_nonsymmv(matrix, eig, eig_vecs, workspace_v);

/* Sort eigenvalues and eigenvectors */
gsl_eigen_nonsymmv_sort(eig, eig_vecs, GSL_EIGEN_SORT_ABS_DESC);

/* print eigenvalues */
printf(" i      sorted eigenvalues     \n");
for(i = 0; i < DIM; i++)
{
c_element = gsl_vector_complex_get(eig, i);
printf("%3d %25.17e + %25.17e * i\n", i, GSL_REAL(c_element), GSL_IMAG(c_element));
}
printf("\n");

/* print eigenvectors */
printf("sorted eigenvectors : [v1 v2 ... vn]\n");
for(i = 0; i < DIM; i++)
{
c_element = gsl_vector_complex_get(eig, i);
printf("eig%2d = %5.1e + %5.1e * i ", i, GSL_REAL(c_element), GSL_IMAG(c_element));
}
printf("\n");
for(i = 0; i < DIM; i++)
{
for(j = 0; j < DIM; j++)
{
c_element = gsl_matrix_complex_get(eig_vecs, i, j);
printf("%15.8e + %15.8e * i, ", GSL_REAL(c_element), GSL_IMAG(c_element));
}
printf("\n");
}
printf("\n");

/* A * [v1 v2 ... vn] - diag[eig1 eig2 ... eign] * [v1 v2 ... vn] == 0 ? */
c_matrix = gsl_matrix_complex_alloc(DIM, DIM);
tmp_matrix = gsl_matrix_complex_alloc(DIM, DIM);
gsl_matrix_complex_set_identity(tmp_matrix);
for(i = 0; i < DIM; i++)
{
for(j = 0; j < DIM; j++)
{
gsl_matrix_complex_set(c_matrix, j, i, gsl_vector_complex_get(eig, i));
GSL_REAL(c_element) = gsl_matrix_get(companion_a, i, j);
GSL_IMAG(c_element) = 0.0;
gsl_matrix_complex_set(tmp_matrix, i, j, c_element);
}
}

gsl_matrix_complex_mul_elements(c_matrix, eig_vecs);

/* alpha = 1, beta = -1 */
GSL_REAL(alpha) =  1.0; GSL_IMAG(alpha) = 0.0;
GSL_REAL(beta)  = -1.0; GSL_IMAG(beta)  = 0.0;

/* compute A * Eigenvectors(A) - eigenvalues * Eigenvectors(A) */
gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, alpha, tmp_matrix, eig_vecs, beta, c_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++)
{
c_element = gsl_matrix_complex_get(c_matrix, i, j);
printf("%g ", GSL_MAX(fabs(GSL_REAL(c_element)), fabs(GSL_IMAG(c_element))));
}
printf("\n");
}
printf("\n");

/* free workspace */
gsl_matrix_complex_free(eig_vecs);
gsl_matrix_complex_free(c_matrix);
gsl_matrix_complex_free(tmp_matrix);
gsl_eigen_nonsymmv_free(workspace_v);

/* free a, x, b */
gsl_matrix_free(companion_a);
gsl_matrix_free(matrix);
gsl_vector_complex_free(eig_vec);
gsl_vector_complex_free(eig);

return 0;
}
```

### Example execution

```\$ ./eigen_nonsymmetric
Companion Matrix(DIM = 5)
0: -0.833333 -0.666667 -0.5 -0.333333 -0.166667
1: 1 0 0 0 0
2: 0 1 0 0 0
3: 0 0 1 0 0
4: 0 0 0 1 0

i       eigenvalues
0  -3.75695199225259346e-01 +   5.70175161011412079e-01 * i
1  -3.75695199225259346e-01 +  -5.70175161011412079e-01 * i
2  -6.70332047603097281e-01 +   0.00000000000000000e+00 * i
3   2.94194556360141191e-01 +   6.68367097443301916e-01 * i
4   2.94194556360141191e-01 +  -6.68367097443301916e-01 * i

i      sorted eigenvalues
0   2.94194556360141690e-01 +   6.68367097443300695e-01 * i
1   2.94194556360141690e-01 +  -6.68367097443300695e-01 * i
2  -3.75695199225259735e-01 +   5.70175161011412523e-01 * i
3  -3.75695199225259735e-01 +  -5.70175161011412523e-01 * i
4  -6.70332047603096504e-01 +   0.00000000000000000e+00 * i

sorted eigenvectors : [v1 v2 ... vn]
eig 0 = 2.9e-01 + 6.7e-01 * i eig 1 = 2.9e-01 + -6.7e-01 * i eig 2 = -3.8e-01 + 5.7e-01 * i eig 3 = -3.8e-01 + -5.7e-01 * i eig 4 = -6.7e-01 + 0.0e+00 * i
1.75051208e-01 +  9.38169222e-02 * i,  1.75051208e-01 + -9.38169222e-02 * i, -9.41339102e-02 + -1.30117121e-01 * i, -9.41339102e-02 +  1.30117121e-01 * i, -1.51221831e-01 +  0.00000000e+00 * i,
2.14158540e-01 + -1.67642800e-01 * i,  2.14158540e-01 +  1.67642800e-01 * i, -8.32690088e-02 +  2.19963418e-01 * i, -8.32690088e-02 + -2.19963418e-01 * i,  2.25592424e-01 +  0.00000000e+00 * i,
-9.19667594e-02 + -3.60901457e-01 * i, -9.19667594e-02 +  3.60901457e-01 * i,  3.36091341e-01 + -7.54134824e-02 * i,  3.36091341e-01 +  7.54134824e-02 * i, -3.36538324e-01 +  0.00000000e+00 * i,
-5.03072154e-01 + -8.38376547e-02 * i, -5.03072154e-01 +  8.38376547e-02 * i, -3.63041358e-01 + -3.50240521e-01 * i, -3.63041358e-01 +  3.50240521e-01 * i,  5.02047194e-01 +  0.00000000e+00 * i,
-3.82615424e-01 +  5.84272896e-01 * i, -3.82615424e-01 + -5.84272896e-01 * i, -1.35776970e-01 +  7.26184061e-01 * i, -1.35776970e-01 + -7.26184061e-01 * i, -7.48952994e-01 +  0.00000000e+00 * i,

A * [v1 v2 ... vn] - diag[eig1 eig2 ... eign] * [v1 v2 ... vn]
0: 2.77556e-16 2.77556e-16 2.39392e-16 2.39392e-16 6.93889e-17
1: 3.05311e-16 3.05311e-16 2.22045e-16 2.22045e-16 2.77556e-17
2: 4.71845e-16 4.71845e-16 2.498e-16 2.498e-16 5.55112e-17
3: 7.77156e-16 7.77156e-16 7.63278e-16 7.63278e-16 3.33067e-16
4: 1.09635e-15 1.09635e-15 6.10623e-16 6.10623e-16 6.66134e-16
```

<- Go back