//#include <config.h>
#include <stdlib.h>
#include <string.h>

#include <gsl/gsl_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_eigen.h>

static void init (){
	  gsl_ieee_env_setup ();
}//init

static void check (unsigned n, double * mat){
	double d = 0.0;
	unsigned n2 = n*n;
	for (unsigned i = 0; i < n2; i++){
		double el = mat[i];
		d += el * el;
	}//i
	printf ("Hilbert-Schmidt-Norm = %f\n",sqrt(d));
}//check

static void eigen (unsigned n, double * mat){
	gsl_vector * eval2 = gsl_vector_alloc(n);
	gsl_eigen_symm_workspace * w1 = gsl_eigen_symm_alloc (n);
	gsl_matrix_view matView = gsl_matrix_view_array (mat, n, n);

	gsl_eigen_symm(&matView.matrix, eval2, w1);

	printf ("Eigenwerte :\n");
	for (unsigned i = 0; i < n; i++){
		printf ("%u\t%f\n",i,gsl_vector_get (eval2, i));
	}//i
}//eigen 

int main (int argc, char * argv[]){
	int rc = 0;

	if (argc > 1){
		FILE * fp = fopen (argv[1],"rb");
		if (fp != NULL){
			unsigned n = 0;

			if (fread(&n,sizeof(n),1,fp) != 1){
				perror("");
				fprintf (stderr,"Kann Datei %s nicht lesen\n", argv[1]);
				rc = 102;
			}
			else {
				printf ("n = %u\n",n);

				double * mat = new double[n*n];
				if (fread(mat,sizeof(mat[0])*n*n,1,fp) != 1){
					perror("");
					fprintf (stderr,"Kann Datei %s nicht lesen\n", argv[1]);
					rc = 103;
				}
				else {	
					check(n,mat);
			
					init();
					//885 statt n : auch loop
					//884 statt n : kein loop
					eigen(n,mat);
				}
				delete [] mat;
			}
			fclose(fp);

		}
		else {
			perror("");
			fprintf (stderr,"Kann Datei %s nicht öffnen\n", argv[1]);
			rc = 101;
		}
	}
	else {
		fprintf (stderr,"Keine Datei angegeben !\n");
		rc = 100;
	}

	printf ("Ende eigenvals, rc = %d\n",rc);
	return rc;
}