/*
/  Author: Haytham Alsaghayer
/  Class: MA 5903
/  Date: Aug 07, 2009
/  Assignment: Conjugate Gradient
*/ 

#include "Header.h"

int main(int argc, char *argv[])
{
	int i, j;				//used in for loops
	//to hold mat file name from stdin
	char matFileName[MAX_FILE_LEN];
	//to hold vec file name from stdin
	char vecFileName[MAX_FILE_LEN];
	
	//from the main problem:  A x = b
	Mat * A;
	Vec * b;
	Vec * x;
	
	gsl_vector_float * gb;		//gsl vector
	gsl_matrix_float * gA;		//gsl matrix
	
	//matrix and vector file streams
	FILE * matInFile;
	FILE * vecInFile;
	
	//get file names from stdin
	scanf("%s%s", matFileName, vecFileName);
	
	//open mat file
	matInFile = fopen(matFileName, "r");
	//open vec file
	vecInFile = fopen(vecFileName, "r");
	
	//read mat from file
	A = readMat(matInFile);
	//read vec from file
	b = readVec(vecInFile);
	
	//check valdity of matrix and vector
	//if (A is not a square matrix) OR (matrix and vector do not conform)
	if((A->n != A->m) || (A->m != b->n))
	{
		printf("ERROR: Matrix is not square OR matrix and vector do not conform\n");
		printf("Program Terminated\n");
		exit(1);
	}

	//allocate gsl matrix
	gA = gsl_matrix_float_alloc(A->n, A->m);
	//copy values from A to gA
	for(i =0; i < A->n; i++)
	{
		for(j =0; j < A->m; j++)
		{
			gsl_matrix_float_set(gA, i, j, A->p[i*A->m + j]);
		}
	}

	//allocate gsl vector
	gb = gsl_vector_float_alloc(b->n);
	//copy values from b to gb
	for(i =0; i < b->n; i++)
	{
		gsl_vector_float_set(gb, i, b->p[i]);
	}

	// start timer
	struct timeval tim;
	gettimeofday(&tim, NULL);
	double t1=tim.tv_sec+(tim.tv_usec/MILLISECOND);
	
	//call Conjugate Gradient function
	x = CG(gA, gb, TOLERANCE, MAX_ITRN);
	
	// end timer
	gettimeofday(&tim, NULL);
	double t2=tim.tv_sec+(tim.tv_usec/MILLISECOND);
	
	//print x
	writeVec(x, stdout);
	//print duration
	printf("\nDuration: %g seconds\n", (t2-t1));
	
	//close streams
	fclose(matInFile);
	fclose(vecInFile);
	
	//free memory resources
	gsl_matrix_float_free(gA);
	gsl_vector_float_free(gb);
	destroyMat(A);
	destroyVec(b);
	destroyVec(x);
	
	return 0;
}