/*
/  Author: Haytham Alsaghayer
/  Class: MA 5903
/  Date: Aug 12, 2009
/  Assignment: QR Eigenvalues
*/ 

#include "Header.h"

int main(int argc, char *argv[])
{
	int i, iLimit;
	Mat * myA;
	Vec * v;
	gsl_eigen_symm_workspace * wSpace;	//work space for eigenvalues
	gsl_matrix * A;
	gsl_vector * eval;
	
	//allocation
	myA	= readMat(stdin);			//read from stdin
	v	= createVec(myA->n);			//allocate vec v
	wSpace	= gsl_eigen_symm_alloc(myA->n);		//allocate work space 
	A	= gsl_matrix_alloc(myA->n, myA->m);	//allocate gsl-matrix
	eval	= gsl_vector_alloc(myA->n);		//allocate gsl-vec
	
	iLimit = A->size1 * A->size2;
	//copy values from myA to A
	for(i =0; i < iLimit; i++)
	{
		A->data[i] = (double) myA->p[i];
	}
	
	//compute eigenvalues using my function
	QREvals(myA, v, MAX_ITRN, TOLERANCE);
	
	//compute eigenvalues using gsl built-in function
	gsl_eigen_symm(A, eval, wSpace);
	//sort eval engenvalues
	qsort (eval->data, eval->size, sizeof(double), gslVecCompare);
	
	//print eigenvalues
	//printf("My function's eigenvalues printed in our Vector format: \n");
	//writeVec(v, stdout);
	//printf("Built in gsl function's eigenvalues: \n");
	//gsl_vector_fprintf(stdout, eval, "%g");
	
	//print
	//my QR		built-in gsl-QR		error
	printf("%-14s   %-14s   %-14s\n",
		"my QR", "gsl-QR", "error");
	for(i =0; i < v->n; i++)
	{
		printf("%-14g   %-14g   %-14g\n",
			v->p[i], eval->data[i],
			max(eval->data[i], (double)v->p[i]) - min(eval->data[i], (double)v->p[i]));
	}
	
	//deallocate (free memory resources)
	destroyMat(myA);
	destroyVec(v);
	gsl_eigen_symm_free(wSpace);
	
	return 0;
}