#include "header.h"

int checksys(gsl_matrix *A)
{
	int i,j;
	int n=A->size1;
	float tmp1;
	float tmp2;
	printf("n====%d\n",n);
	for(i=1;i<n;i++)
	{
	   for(j=i+1;j<n;j++)
	   {
		tmp1=gsl_matrix_get(A,i,j);
		tmp2=gsl_matrix_get(A,j,i);
		if(tmp1!=tmp2) return 0;
	   }
	}
	return 1;
}

void VectorSub(gsl_vector *a,gsl_vector *b,gsl_vector *ans)
{
	gsl_vector_memcpy(ans,a);
	gsl_vector_sub(ans,b);
}

int QREvals(gsl_matrix **A,gsl_vector **v, int maxn, float Tol)
{
	int i=0,j=0,k=0;
	int n;
	double tmp;
	double err;
	int count=0;
	gsl_matrix *A0;
	gsl_matrix *Q;
	gsl_vector *v0;
	gsl_vector *r;
	n=(*A)->size1;
	Q=gsl_matrix_alloc (n,n);
	A0=gsl_matrix_alloc (n,n);
	v0=gsl_vector_alloc (n);
	r=gsl_vector_alloc (n);
///////////check symmetric 
	if (!checksys(*A))
	{
		printf("the input matrix is not symmetric!\n");
		printf("program stopped!\n");
		return 0;
	}
///////////A is symmetric
//////////QR algorithm begin
for(count=0;count<maxn;count++)
{
	gsl_matrix_memcpy(A0,*A);
	gsl_vector_memcpy(v0,*v);
/////////////////////////////QR decomp///////////////////////
	gsl_linalg_QR_decomp(*A,*v); // this function return a wrong answer comparing to using R to check.
/////////////////////////////calculate err///////////////////////
	VectorSub(*v,v0,r);         // r=v-v0;
	if(gsl_blas_dnrm2(r)<Tol)
		break;
/////////////////////////////Get Q//////////////////////
	for(j=0;j<n;j++)
	{
		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  
		}
	}
 /* Compute A0 = t(Q)*A */
       gsl_blas_dgemm (CblasTrans, CblasNoTrans,
                       1.0, Q, *A,
                       0.0, A0);
       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;        // error
	int maxite;       //maximum iterate
	int n,m;            // dimentions of A;
	int flag=0;
/////////////////// input begins
	//////////// read in the matrix from the standard input stream.
	scanf("%d",&n);
	scanf("%d",&m);
	A=gsl_matrix_alloc(n,m);   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	B=gsl_matrix_alloc(n,m);   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	gsl_matrix_fscanf (stdin,A);
	gsl_matrix_memcpy(B,A);
/////////////////// call the CG function to solve the equation
        x=gsl_vector_alloc(n);  //!!!!!!!!!!!!!!!!!!!!
	gsl_vector_set_zero(x); //x=0;
	tol=TOL;                      //control the error
	maxite=MAXITE;                //control the iterate
	flag=QREvals(&A,&x,maxite,tol);
////////////////////compare with build-in function
        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("eval by using gsl function:\n");
	gsl_vector_fprintf(stdout,eval,"%lf");
/////////////////output begins
	if(flag==1)
	{
		printf("QR Eigenvalues of Matrix A is: \n");
		gsl_vector_fprintf(stdout,eval,"%lf");
	}
////////////////
	gsl_vector_free(x);
	gsl_matrix_free(A);
	return 0;
}