#define RND ((double)rand()/RAND_MAX)
#include <stdio.h>
#include <time.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#define TIME ((double)clock())/CLOCKS_PER_SEC

int main(int argc, char** argv){
	printf("argc=%d\n",argc);
	int n=(argc>1? atoi(argv[1]):10);
	printf("n=%d\n",n);
	
//	gsl_matrix *A = gsl_matrix_alloc(n,n);
	gsl_matrix *B = gsl_matrix_alloc(n,n);
//	gsl_matrix *Q = gsl_matrix_alloc(n,n);
//	gsl_matrix *R = gsl_matrix_alloc(n,n);
	gsl_matrix *V = gsl_matrix_alloc(n,n);
	gsl_vector *e = gsl_vector_alloc(n);
	gsl_vector *t = gsl_vector_alloc(n);
	gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc(n);

	for(int i=0;i<n;i++) for(int j=i;j<n;j++){
//		double x= (i==j? (double)(i+1):0);
//		gsl_matrix_set(A,i,j,x);
//		gsl_matrix_set(A,j,i,x);
		double rnd=RND;
		gsl_matrix_set(B,i,j,rnd);
		gsl_matrix_set(B,j,i,rnd);
	}

//	gsl_linalg_QR_decomp(B,t);
//	gsl_linalg_QR_unpack(B,t,Q,R);

//	gsl_blas_dgemm(CblasTrans,CblasNoTrans,1.0,Q,A,0.0,V);
//	gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,V,Q,0.0,A);

//	double time=TIME;
//	gsl_eigen_symmv(A,e,V,w);
	gsl_eigen_symmv(B,e,V,w);
//	time=TIME-time;

//	for(int i=0;i<n;i++) printf("e(n)=%g\n",gsl_vector_get(e,i));
//	printf("time\t = %g\n",time);
	printf("e(n)=%g\n",gsl_vector_get(e,n-1));
	printf("e(1)=%g\n",gsl_vector_get(e,0));

	gsl_eigen_symmv_free(w);
//	gsl_matrix_free(A);
	gsl_matrix_free(V);
	gsl_vector_free(e);
	
	return 0;
}