9 de julho de 2009

Using GSL to calculate Nonsymmetric matrices eigenvalues

Here is a example of using GSL - Gnu Scientific Library - to calculate the eigenvalues of a nonsymmetric matrix.
This is a simple example write in C. GSL uses BLAS and can be optimized with ATLAS, so can be used as a real tool for heavy calculations. To get know more about GSL visit http://www.gnu.org/software/gsl/ or the manual page.

Hint: Beware when using GSL or our own linear algebra routine, remember that greats matrices requires greats responsibilities, so keep track of precision errors. Find precision errors will be a great excuse to know more about multiple precision GMP .

#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
int main (void)
int a=2; //matrix a x a
int i,j;
double cg[a][a];

// Defining matrix
cg[0][0] = 2.;
cg[0][1] = 3.;
cg[1][0] = -1.;
cg[1][1] = 1.;

gsl_matrix *TOTAL = gsl_matrix_alloc (a,a); //Alloc GSL matrix TOTAL

//Setting TOTAL
for (i = 0; i < a; i++)
for (j = 0; j < a; j++)
gsl_matrix_set (TOTAL, i, j, cg[i][j]);

gsl_vector_complex *eval = gsl_vector_complex_alloc (a); //Vector whose entries are the eigenvalues
gsl_eigen_nonsymm_workspace *K = gsl_eigen_nonsymm_alloc (a); //alloc work space for GSL nonsymmetric calculations
gsl_eigen_nonsymm (TOTAL,eval,K);
gsl_eigen_nonsymm_free (K);
gsl_matrix_free (TOTAL);

for (i = 0; i < a; i++) /* OUT OF RANGE ERROR */
gsl_complex eval_i = gsl_vector_complex_get (eval, i);
printf ("%g + %g*i\n",GSL_REAL(eval_i),GSL_IMAG(eval_i) );
return 0;