#include <iostream.h>
#include <vector.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>

//the function that solves linear functions and print the results
void solve(vector<double>, vector<double>);

int main(){
	vector<double> x, matrix_elements;
	//The vector that contains diagonal and off-diagnonal elements
	matrix_elements.push_back(2);
	matrix_elements.push_back(4);
	//The vector that contains x vector
	x.push_back(1);
	x.push_back(2);
	x.push_back(3);
	//comparison of the outputs from two methods
	solve(x, matrix_elements);
	return 0;
}
//comparison of the outputs from two methods
void solve(vector<double> x_value, vector<double> matrix_value){
	//The dimension of the matrix 
	int length_diag = x_value.size();
	//Two copies of the matrices
	gsl_matrix *varMatrix, *copy_varMatrix;
	//Two copies of x vectors
	gsl_vector *xV = gsl_vector_calloc(length_diag), *copy_xV = gsl_vector_calloc(length_diag);
	for(int i=0;i<length_diag;i++){
		gsl_vector_set(xV,i,x_value[i]);
		gsl_vector_set(copy_xV, i, x_value[i]);
	}
	varMatrix = gsl_matrix_alloc(length_diag, length_diag);
	copy_varMatrix = gsl_matrix_alloc(length_diag, length_diag);
	cout<<"********Comparison 1: The calculation for a symetric matrix************"<<endl;
	//set up symetric matrix
	for(int i=0;i<length_diag;i++)
		 for(int j=0;j<length_diag;j++){
			if(i==j){
				gsl_matrix_set(varMatrix,i,j,matrix_value[0]+matrix_value[1]);
				gsl_matrix_set(copy_varMatrix,i,j,matrix_value[0] + matrix_value[1]);
			}
			else {
				gsl_matrix_set(varMatrix,i,j,matrix_value[1]);
				gsl_matrix_set(copy_varMatrix,i,j,matrix_value[1]);
			}
		}
	cout<<"The input symetric matrix is "<<endl;
	for(int i=0;i<length_diag;i++){
		for(int j=0;j<length_diag;j++){
			cout<<gsl_matrix_get(varMatrix,i,j)<<'\t';
		}
		 cout<<endl;
	}							
	gsl_linalg_cholesky_decomp(varMatrix);
    cout<<"The matrix after cholesky decomposition is"<<endl;
	for(int i=0;i<length_diag;i++){
		for(int j=0;j<length_diag;j++){
			cout<<gsl_matrix_get(varMatrix,i,j)<<'\t';
		}
		cout<<endl;
	}
		 
	gsl_linalg_cholesky_svx(varMatrix,xV);
	cout<<"Comparison of the outputs from gsl_linalg_cholesky_svx and gsl_linalg_HH_svx"<<endl;
	cout<<"The output from gsl_linalg_cholesky_svx"<<endl;
	for(int i=0;i<length_diag;i++){
		cout<<gsl_vector_get(xV,i)<<'\t';
	}
	cout<<endl;
	cout<<"The output from gsl_linalg_HH_svx"<<endl;
	gsl_linalg_HH_svx(copy_varMatrix,copy_xV);
	for(int i=0;i<length_diag;i++){
		cout<<gsl_vector_get(copy_xV,i)<<'\t';
	} cout<<endl;
	cout<<"\n\n\n"<<endl;
	cout<<"********Comparison 2: The calculation for a matrix with upper triangle elements being zeroes************"<<endl;
	// reset x vectors
	for(int i=0;i<length_diag;i++){
		gsl_vector_set(xV,i,x_value[i]);
		gsl_vector_set(copy_xV, i, x_value[i]);
	}
								
	//set up a matrix with upper triangle elements being zeroes
	for(int i=0;i<length_diag;i++)
		for(int j=0;j<length_diag;j++){
			if(i==j){
				gsl_matrix_set(varMatrix,i,j,matrix_value[0]+matrix_value[1]);
				gsl_matrix_set(copy_varMatrix,i,j,matrix_value[0] + matrix_value[1]);
			}
			else if(i>j){
				gsl_matrix_set(varMatrix,i,j,matrix_value[1]);
				gsl_matrix_set(copy_varMatrix,i,j,matrix_value[1]);
			}
		
			//set upper triangle to be zeroes
			else {
				gsl_matrix_set(varMatrix,i,j,0);
				gsl_matrix_set(copy_varMatrix,i,j,0);
			}
		
		}
	cout<<"The input matrix is "<<endl;
	for(int i=0;i<length_diag;i++){
		for(int j=0;j<length_diag;j++){
			cout<<gsl_matrix_get(copy_varMatrix,i,j)<<'\t';
		}
		cout<<endl;
	}
	gsl_linalg_cholesky_decomp(varMatrix);
	cout<<"The matrix after cholesky decomposition is"<<endl;
	for(int i=0;i<length_diag;i++){
		for(int j=0;j<length_diag;j++){
			cout<<gsl_matrix_get(varMatrix,i,j)<<'\t';
		}
		cout<<endl;
	}
	gsl_linalg_cholesky_svx(varMatrix,xV);
	cout<<"Comparison of the outputs from gsl_linalg_cholesky_svx and gsl_linalg_HH_svx"<<endl;
	cout<<"The output from gsl_linalg_cholesky_svx"<<endl;
	for(int i=0;i<length_diag;i++){
		cout<<gsl_vector_get(xV,i)<<'\t';
	}
	cout<<endl;
	cout<<"The output from gsl_linalg_HH_svx"<<endl;
	gsl_linalg_HH_svx(copy_varMatrix,copy_xV);
	for(int i=0;i<length_diag;i++){
		cout<<gsl_vector_get(copy_xV,i)<<'\t';
	}
	cout<<endl;

}