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;
}