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