/**  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>

#include "version.h"

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

	fprintf(stdout, "%s\n", ESPACKVersion);

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