#include "header.h"

int QREvals (gsl_matrix *A,gsl_vector *seigen,int n,float Tol)
{       

///////////////////////////// allocate 
	int i,j,count,max=1000,ss=0;
	double tempa0,tempa1,temp0; //get the diagonal element of A and A0
	double static error;
	gsl_matrix *tempA,*invR,*Q,*A0,*S,*AA;  //S is used to check the Q and R is right
	gsl_vector *tau,*a0,*a1,*eigen;
	tau=gsl_vector_alloc(n);
	tempA=gsl_matrix_alloc(n,n); 

	invR=gsl_matrix_alloc(n,n); 
	A0=gsl_matrix_alloc(n,n); 
	AA=gsl_matrix_alloc(n,n);
	S=gsl_matrix_alloc(n,n); 
	Q=gsl_matrix_alloc(n,n);  
	
	gsl_matrix_set_all(S,0.0);
	a0=gsl_vector_alloc(n); 
	a1=gsl_vector_alloc(n); 
	eigen=gsl_vector_alloc(n); 

	 for(count=0;count<max;count++)	{
	/////////////////////////qr decomposition
		gsl_matrix_memcpy(AA,A);
		gsl_linalg_QR_decomp(AA,tau);
		gsl_matrix_memcpy(tempA,AA); // copy A into tempA
		for (i=0;i<n;i++) {
			for (j=0;j<n;j++) 
			{
				if (j<i)  gsl_matrix_set(tempA,i,j,0.0);   // set up the matrix R in tempA
			}
		}
		gsl_matrix_set_identity(invR);
		gsl_blas_dtrsm(CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit,1.0,tempA,invR); // get the inverse of R in invR.
	
		//gsl_matrix_fprintf(stdout,tempA,"%lf");
		//gsl_matrix_fprintf(stdout,A,"%lf");
		//printf("###################\n");
		gsl_matrix_set_all(Q,0.0);
		gsl_blas_dgemm (CblasNoTrans,CblasNoTrans,1.0,A,invR,1.0,Q); // get the matrix Q;
		gsl_matrix_set_all(A0,0.0);
		gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,tempA,Q,1.0,A0); // get the matrix A0;
		//gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,Q,tempA,1.0,S); 
		//gsl_matrix_fprintf(stdout,A0,"%lf");
	/////////////get the diagonal value in vector a1, a0, the norm 2 error in error	
		for (i=0;i<n;i++) {
			for (j=0;j<=n;j++)
	  		{	
				if (i==j) {
				tempa1=gsl_matrix_get(A,i,j);
				gsl_vector_set(a1,i,tempa1);
				tempa0=gsl_matrix_get(A0,i,j);
				gsl_vector_set(a0,i,tempa0);
				}
			}
		}
		gsl_vector_sub(a1,a0);  // use a1 to store the error
		error=gsl_blas_dnrm2(a1);
		ss++;	
		if (error<Tol) break; 
		gsl_matrix_memcpy(A,A0);
	 }
///////////////////////////////////////////////////////////////////////////////////////////////
	for(j=0;j<n;j++)
		{
			gsl_vector_set(eigen,j,gsl_matrix_get(A,j,j));    //diag
		}
	gsl_vector_fprintf(stdout,eigen,"%lf");
	gsl_sort_vector(eigen);           
	for(i=0;i<n;i++)
	{
		 gsl_vector_set(seigen,i,gsl_vector_get(eigen,n-1-i)); // from largest to smallest
	}
	return ss;
}