#include "header.h"


int QREvals(gsl_matrix **A,gsl_vector **v, int maxn, float Tol)
{
	int i=0,j=0,k=0;
	int n;
	int count=0;
	gsl_matrix *A0;
	gsl_matrix *Q;
	gsl_vector *v0;
	gsl_vector *r;

	n=(*A)->size1;
	Q=gsl_matrix_alloc (n,n);	// initialize the n by n matrix: Q
	A0=gsl_matrix_alloc (n,n);	// initialize the n by n matrix: A0
	v0=gsl_vector_alloc (n);	// initialize the n by 1 vector: v0
	r=gsl_vector_alloc (n);		// initialize the n by 1 vector: r

	if (!check(*A))	//check symmetry 
	{
		printf("The matrix A is not symmetric!\n");
		return 0;
	}

	for(count=0;count<maxn;count++)
	{
		gsl_matrix_memcpy(A0,*A);
		gsl_vector_memcpy(v0,*v);

		gsl_linalg_QR_decomp(*A,*v); // QR

		VectorSub(*v,v0,r);         // r=v-v0
		if(gsl_blas_dnrm2(r)<Tol)
			break;

		for(j=0;j<n;j++)	// find Q
		{
			for(k=0;k<n;k++)
			{
				if(j==k) gsl_matrix_set(Q,j,k,gsl_vector_get(*v,j));    //diag
				if(j<k)  gsl_matrix_set(Q,j,k,0);			// upper trig 
				if(j>k)  gsl_matrix_set(Q,j,k,gsl_matrix_get(*A,j,k));	// lower trig  
			}
		}

	       gsl_blas_dgemm (CblasTrans, CblasNoTrans,1.0, Q, *A,0.0, A0);	// Compute A0 = t(Q)*A 
	       gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,1.0, A0,Q,0.0, *A);
	}
	return 1;
}


int main(void)
{
	gsl_matrix *A;
	gsl_matrix *B;
	gsl_vector *x;
	float tol; 
	int maxite;     //maximum iterate
	int n,m;          // dimentions of A;
	int flag=0;

	scanf("%d",&n);	// read in matrix from the standard input
	scanf("%d",&m);
	A=gsl_matrix_alloc(n,m);   // A is n by m
	B=gsl_matrix_alloc(n,m);   // B is n by m
	gsl_matrix_fscanf (stdin,A);	// read A
	gsl_matrix_memcpy(B,A);	// assign A to B

      x=gsl_vector_alloc(n);   // x is n by 1
	gsl_vector_set_zero(x); // initialize x as a zero vector

	flag=QREvals(&A,&x,MAXITE,TOL);	// build-in function, return an integer
      gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (n);
      gsl_vector *eval = gsl_vector_alloc (n);
      gsl_matrix *evec = gsl_matrix_alloc (n, n);
      gsl_eigen_symmv (B, eval, evec, w);
	printf("e-values by using gsl function:\n");
	gsl_vector_fprintf(stdout,eval,"%lf");

	if(flag==1)	// A is symmetric
	{
		printf("QR e-values of Matrix A is: \n");
		gsl_vector_fprintf(stdout, x,"%lf");
	}

	gsl_vector_free(x);
	gsl_matrix_free(A);
	return 0;
}