#include "Header.h"

//computes the eigenvalues of A, return the ordered (from largest
// to smallest) eigenvalues in v, with a maximum of n iterations,
// untill it satisfies some error estimate less than Tol
//
//A and v must be allocated and with the right size
int QREvals(Mat * myA,  Vec  * v, int n, float Tol)
{
	//check that myA is symmetric (A = A^T)
	//  if myA is not symmetric, exit
	if(!(isMatSymm(myA)))
	{
		printf("Matrix is not symmetric.\n");
		printf("Program aborted.\n");
		exit(0);
	}
	int i, iLimit;
	gsl_vector 	* tau;
	gsl_matrix 	* Q;
	gsl_matrix	* R;
	gsl_matrix	* A;
	
	//allocation
	A	=	gsl_matrix_alloc(myA->n, myA->m);
	tau	=	gsl_vector_alloc(A->size1);
	Q	=	gsl_matrix_alloc(A->size1, A->size1);
	R	=	gsl_matrix_alloc(A->size1, A->size2);
	
	iLimit = A->size1 * A->size2;
	//copy values from myA to A
	for(i =0; i < iLimit; i++)
	{
		A->data[i] = (double) myA->p[i];
	}
	
	//do QR algorithm
	// while (maxItr (n) not reached) OR (A->data[0] is bigger than tolerance (Tol))
	for(i =0; (i < n) && (1 -(A->data[0]-((int)A->data[0])) > Tol); i++)
	{
		//QR decomp
		gsl_linalg_QR_decomp(A, tau);
		//QR unpack to Q and R
		gsl_linalg_QR_unpack(A, tau, Q, R);
		//A = R*Q
		gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, R, Q, 0.0, A);
	}
	
	//copy result to v
	for(i =0; i < A->size1; i++)
	{
		v->p[i] = (float) A->data[i*A->size1+i];
	}
	
	//sort data
	qsort (v->p, v->n, sizeof(float), vecCompare);
	
	//deallocate (free memory resources)
	gsl_vector_free(tau);
	gsl_matrix_free(A);
	gsl_matrix_free(Q);
	gsl_matrix_free(R);
	
	return 1;
}