#include "header.h"

// This function is using QR decompesition to find eigenvalues of a symmetric matrix.
int QREvals(gsl_matrix *A,gsl_vector *e, int m, float Tol)
{
	int i,j,k,c,n;
	double tmp,err;
	gsl_matrix *A0,*Q,*R;
	gsl_vector *e0,*v,*r;
	n=A->size1;
	A0=gsl_matrix_alloc (n,n);	
	Q=gsl_matrix_alloc (n,n);
	R=gsl_matrix_alloc (n,n);
	e0=gsl_vector_alloc (n);
	v=gsl_vector_alloc (n);	
	r=gsl_vector_alloc (n);

	if (!ifsym(A))
	{
		printf("the input matrix is not symmetric!\n");
		printf("That may cause the result way from the ture eigenvalue of matrix.\n");
		printf("Hence stop calculating!\n");
		return 0;
	}//symmetric is important here for QR algorithm

	for(j=0;j<n;j++)
	{
		gsl_vector_set(e0,j,gsl_matrix_get(A,j,j));
	}
	//////////QR algorithm begin
	for(c=0;c<m;c++)
	{
		gsl_matrix_memcpy(A0,A);
	
		gsl_linalg_QR_decomp(A,v); //QR decomp using the gsl function
		
		//next loop is to get the matrix R(which is a upper trig with diag not all zero).
		for(j=0;j<n;j++)
		{
			for(k=0;k<n;k++)
			{
				if(j>k) gsl_matrix_set(R,j,k,0); //lower trig
				else    gsl_matrix_set(R,j,k,gsl_matrix_get(A,j,k));    //diag & upper trig

			}
		}
//Then get Q(orthognal matrix)
	gsl_matrix_memcpy(A,A0);
	gsl_blas_dtrsm (CblasRight,CblasUpper,CblasNoTrans,CblasNonUnit,1.0,R,A);
	gsl_matrix_memcpy(Q,A);	
	
//Followed is to find next A(called A') and set to A
	 
	       gsl_matrix_memcpy(A,A0);
	       gsl_blas_dgemm (CblasTrans, CblasNoTrans,
		               1.0, Q, A,
		               0.0, A0);
	 
	       gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
		               1.0, A0,Q,
		               0.0, A);
	//Get eigenvalues(which is a vector) for A,by using the algorithm in readme.
		for(j=0;j<n;j++)
		{
			gsl_vector_set(e,j,gsl_matrix_get(A,j,j));    //diag(as mentioned in readme, the diag of A
		}
		//calculate err///////////////////////
		vsub(e,e0,r);  //e0 the previous one and e recent one       
		if(gsl_blas_dnrm2(r)<Tol)
			break;
		gsl_vector_memcpy(e0,e);
	}
	return c;
}