#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);
		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);
	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));

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