Last Update: Aug. 29, 2008
/***********************************************/
/* 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;
}
$ ./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