source: util/src/diagonalize.c@ a1c448

/**  Diagonlizes a 3x3 matrix.
 * Most of here comes from the GSL example
 */

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

int main(int argc, char **argv)
{
        double matrix[9];
        int i;

        if (argc < 9) {
                printf("Usage: %s <9 entries of 3x3 matrix>\n", argv[0]);
                return 1;
        } else {
                for(i=0;i<9;i++)
                        matrix[i] = atof(argv[i+1]);
        }

        gsl_matrix_view m 
                = gsl_matrix_view_array (matrix, 3, 3);
     
        gsl_vector *eval = gsl_vector_alloc (3);
        gsl_matrix *evec = gsl_matrix_alloc (3, 3);
     
        gsl_eigen_symmv_workspace * w = 
          gsl_eigen_symmv_alloc (4);
        
        gsl_eigen_symmv (&m.matrix, eval, evec, w);
     
        gsl_eigen_symmv_free (w);
     
        gsl_eigen_symmv_sort (eval, evec, 
                                 GSL_EIGEN_SORT_ABS_ASC);
        
        {
          for (i = 0; i < 3; i++)
            {
              double eval_i 
                  = gsl_vector_get (eval, i);
              gsl_vector_view evec_i 
                  = gsl_matrix_column (evec, i);
     
              fprintf (stderr, "eigenvalue = %g\n", eval_i);
              fprintf (stderr, "eigenvector = \n");
              gsl_vector_fprintf (stderr, 
                                     &evec_i.vector, "%g");
            }
        }
        fprintf (stdout, "eigenvalues: %lg\t%lg\t%lg\n", gsl_vector_get (eval, 0), gsl_vector_get (eval, 1), gsl_vector_get (eval, 2));
     
        gsl_vector_free (eval);
        gsl_matrix_free (evec);
     

        return 0;
}